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1.   INTRODUCTION 

Unsteady  flow  in  the  passages  of  wave  rotor  devices  can  adequately  be 
modelled  on  a  one-dimensional  basis.   However,  this  modelling  can  be  quite 
involved  due  to  the  peculiar  characteristics  typical  of  wave  rotor  type  flows. 
The  numerical  calculation  has  to  provide  approximate  solutions  of 
time-dependent  compressible  fluid  flow  problems  which  involve  discontinuities 
and  strong  wave  interactions.   Ref .  (1)  lists  three  criteria  which  such 
approximate  solutions  should  satisfy  simultaneously:   (i)  the  solution  must  be 
reasonably  accurate  in  smooth  regions  of  the  flow.   Continuous  waves 
(rarefaction  waves,  compression  waves)  should  propagate  at  the  correct  speed 
and  should  maintain  the  correct  shape  which  involves  steepening  or  spreading 
at  the  correct  rate;   (ii)  discontinuities  which  are  transported  along 
characteristics  (gradient  discontinuities,  contact  surfaces),  should  be 
modelled  by  sharp  and  discrete  jumps,  and  should  be  transported  at  the  correct 
speed;  and   (iii)  nonlinear  discontinuities  such  as  shocks  should  be  computed 
stably  and  accurately. 

In  addition,  the  complex  pattern  of  shock  waves  and  contact  surfaces  that 
could  evolve  in  wave  rotor  devices  precludes  the  use  of  numerical  methods 
which  rely  on  either  some  typt  of  artificial  viscosity  or  a  special  tri-itment 
of  discontinuities.   Such  methods  would  quickly  become  quite  impractic.il  for 
this  application  due  to  programming  difficulties  and  cost  of  execution. 

Computation  of  such  solutions  has  generally  been  carried  out  by  solving  a 
set  of  finite  difference  equations  which  approximate  the  guverning 
differential  equations  of  flow.   All  such  schemes  inherently  have  a  finite 
amount  of  dissipation  as  well  as  dispersion  of  the  wave  modes  they  model,  and 
it  is  difficult  to  construct  difference  schemes  which  simultaneously  satisfy 
the  criteria  given  above.   Stability  problems  may  also  be  an  added  concern  for 
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these  schemes. 

In  view  of  the  foregoing,  an  alternative  approach  to  solving  wave  rotor 
type  flows  was  sought,  and  the  purpose  of  this  report  is  to  describe  such  a 
scheme  along  with  some  results.   The  scheme  is  known  variously  as  Glimm's 
method,  the  Random  Choice  Method  (RCM)  or  the  piecewise  sampling  method.   The 
method  evolved  from  a  constructive  proof  of  the  existence  of  solutions  to 
systems  of  nonlinear  hyperbolic  conservation  laws  given  by  Glimra  (Ref.  2). 
Chorin  (Refs.  3  and  4)  developed  the  scheme  into  an  effective  numerical  tool 
for  gas  dynamic  applications,  with  emphasis  on  detonation  combustion  problems 
and  reacting  gas  flows.   Although  the  RCM  computes  solutions  on  a  fixed  grid, 
it  is  not  a  difference  scheme,  utilizing  solutions  of  locally  defined  Riemann 
problems  as  the  basic  building  blocks  for  the  global  solution.   Each  of  the 
local  Riemann  problems  (defined  in  more  detail  in  section  2)  provides  an 
analytically  exact  elementary  similarity  solution.   By  means  of  a  suitable 
sampling  procedure,  usually  of  a  pseudo-random  or  quasi-random  nature,  the 
similarity  solutions  are  superposed  to  construct  the  approximate  solution  to 
the  equations. 

With  an  appropriate  sampling  technique,  the  RCM  in  one  dimension  is 
possibly  superior  to  any  finite  difference  scheme  in  meeting  the  criteria 
established  above. 


2.   METHOD 


2.1.   Solution  Procedure 


The   method   models    the    one-dimensional,    compressible,    inviscid    Euler 

equations,  expressed  in  conservation  form  as 

3U    3F(U)    n  , 

-*—  +  —rr — -  -  0   ,   where 
ot      dx 


U(x,t)  = 


pu 
E 


and     F(U)  = 


pu-  +  p 
(E  +  p)„ 


(1) 


Here   E   is  the  total  energy  per  unit  volume  and  may  be  expressed  as  (tor  a 
polytropic  gas) 

E  =  p  e  +  —  ,ju2   ,   e  A  internal  energy  per  unit  mass 

p   is  the  density,   p   is  pressure  and   u   is  velocity  in  the  one  space 
dimension  being  considered  here.   With  initial  data  specified  in  the  form 

U(x,0)  =  9(x)  , 
an  initial  value  problem  is  defined  for  the  Euler  equations.   The  simplest 
initial  value  problem  for  which  discontinuities  appear  is  the  Riemann  probl 
to  find  the  gas  flow  resulting  from  an  initial  state  in  which  the  gas  on  the 
right  of  an  'origin'  is  in  a  constant  state,  and  the  gas  on  the  left  is  in 
another  constant  state,  i.e., 

UL   ,   x  <  0 


era 


•Kx)  = 


with 


UR   ,   x  >  0 


U 


L,R  = 


^L,R 

(^U)L,R 
EL,R 


where  the  subsripts   L  and   R  denote  the  left  and  right  sides  of  the 
'origin',  here  arbitrarily  prescibed  at   0   .   That  is,  the  Riemann  problem 
consists  of  prescribing  constant  initial  data  on  either  side  of  an  origin 
where  a  jump  discontinuity  exists.   As  mentioned  before,  the  solution  of  the 
problem  constitutes  a  basic  building  block  of  the  random  choice  method.   A 
special  case  of  the  Rieraann  problem  in  which  u^,  =  ur  =  0   is  often  referred 
to  as  the  shock  tube  problem.   The  answer  to  the  problem  is  that  there  are 
four  possible  types  of  subsequent  flow,  depending  on  the  inequalities  in  the 
left  and  right  side  data  prescribed.   Thus,  in  both  directions  from  the 
origin,  a  shock  or  a  centered  rarefaction  wave  may  propagate,  giving  rise  to 
the  above  mentioned  four  different  possibilities.   Fig.  (1)  illustrates  the 
special  case  of  shock  tube  type  flow  and  the  evolution  of  the  wave  pattern. 
Fig.  (2)  shows  the  simple  fixed  Cartesian  grid  set  up  for  the  method. 
Let   /ax   be  a  spatial  increment  and   At   a  time  increment.   The  solution  is  to 
be  evaluated  at  time   (n  +  1 ) At   ,   n   being  a  non-negative  integer,  at 
spatial  increments   i Ax   ,   i  =  1,2,3,  .  .  .   The  initial  data  is  prescribed 
for  each  time  step  at   nAt   in  a  piecewise  constant  manner  i.e.,  it  consists 
of  intervals  of  length   Ax  where  the  data  is  constant,  separated  by  jump 
discontinuities : 

U(x  ,  nAt)  =  U\*   ,  (i— )Ax  <  x  <  (i+^-)Ax 

The  solution  at  time   (n+l)At   then  is  required  to  have  the  same  property, 
i.e.,  it  is  piecewise  constant  over  an  interval   Ax   ,  and  it  serves  as  the 
initial  data  for  the  next  time  step: 

U(x,(n+l)At)  =  Uj+1  ,   (i  --j)Ax  <  x  <  (i-^-)Ax 


This  procedure  defines  a  sequence  of  local  Riemann  problems  to  be  solved  at 
each  time  level.   On  the  grid  shown  in  Fig.  2,  for  example,  initial  data  would 
be  specified  at  points  L,  3,  5  .....  setting  up  a  succession  of  Riemann 
problems  defined  by  each  pair  of  states  (1,3)  ,  (3,5)  ,  (5,7),  with  the 
discontinuities  at  the  midpoint  of  each,  i.e.,  at  2,  4,  6,  ....  etc.   If  the 
time  step  increment   At   is  calculated  such  that 

At  <  o.(Ax).  max  ( I u . I +a . )   ,  with 

0  <  o  <  i 

then  the  waves  generated  at  the  discontinuities  of  adjacent  Riemann  problems 
will  not  interact,  as  shown  schematically  in  Fig.  2. 

Each  of  the  local  Riemann  problems  yields  an  exact  analytical  solution, 
with  the  resulting  wave  structure  a  particular  combination/variation  of  the 
general  structure  shown  in  Fig.  3. 

In  the   x-t   plane,  the  solution  to  a  Riemann  problem  consists  of 
essentially  four  regions  connected  by  three  waves.   Thus  states  I  and  TV  are 
the  prescribed  left  and  right  states  for  the  problem,  and  states  II  and  III 

are  the  'starred'  middle  states  separated  by  a  slip  line  or  contact 

dx 
discontinuity  —  =  u*  .   The  velocity,   u   ,  and  pressure,   p   ,  are 

continuous  across  the  contact,  but  ,j       in  general  is  not.   Thus   u^*  =  ur*   , 

PL*  =  PR*   a°d   PL*  *   PR*   •   sl'b   »   s2»b   and   sl»f   »  s2»f   represent 
respectively  the  backward  and  forward  facing  waves  generated  at  the  point  of 
discontinuity  and  may  be  either  shocks  or  rarefaction  waves. 

Still  referring  to  Fig.  3,  it  is  seen  that  at  a  time   nAt  <  t  <  (n+l)At, 
the  exact  solution  of  the  local  Riemann  problem  for  the  interval   [(i-l)Ax   , 
iAx]   may  actually  consist  of  several  distinct  states.   Consider  now  a 


translation  of  each  interval   [(i-l)ax,  iAx]   to 


ax 


T  '   "7 


5] 


such  that 


the  discontinuity  (i.e.,  the  point  from  which  the  waves  are  generated)  is 

centered  at  a  zero  origin.   Let   0  be  the  value  of  a  random  variable,  defined 
over  the  interval   [—=-  ,  +  -=-]   ,  and  let 


Ax  .  _  .    Ax 
C  =  UAx   ,   i.e. 2        ~2 


,n+ 


Also,  define   U"      (x,t)   ,   nAt  <  t  <  (n+l)At  ,  to  be  the  exact  solution  to 
exact 

each  Riemann  problem.   Using  the  value  of  c,      to  fix  a  point  in  the  interval 
ax  of  each  Riemann  problem,  the  exact  solution  at  that  point  is  determined 
and  assigned  Lo  either  the  left  or  the  right  grid  point,  depending  on  whether 
£  is   <  or   >   0   .   Thus,  if  the  point  fixed  by  £,      is   P*   (Fig.  3),  the 
exact  solution  to  the  Riemann  problem  at  that  sampled  location  is  assigned  to 
the  grid  point  on  the  right  and  if  the  sampled  point  is   P"   ,  the  solution  at 
that  location  is  assigned  to  the  grid  point  on  the  left,  i.e.,  for  a  typical 

interval   [(i-l)ax,   iax]   , 

if    C  <  0   ,   U\+[   =   Un+      (C  ,  t) 
l-l       exact 

j  j  c         -   v  n  ,,ri+l     1Tn+     ,  _     . 

and  if    %   >  0   ,   U.    =U       (t,,t) 

l       exact 

It  is  seen  immediately  that  although  the  solutions  are  computed  on  a 

grid  in  this  method,  it  is  not  a  differencing  scheme.   Also,  instead  of  using 

a  weighted  average  of  the  Riemann  problem  solution  to  arrive  at  the  solution 

for  a  grid  pointt,  the  RCM  samples  a  particular  value  from  an  explicit  wave 


t   The  Godunov  method,  for  example  implements 


n+1    1   /*(i+2>Ax 


n+1   j_  n1+2 

1    '  ^J    (l-L 


)Ax 


Un+     (x,t)dx 
exact 


solution,  thus  eliminating  the  smoothing  out  of  wave  transport  and  interaction 
information  inherent  in  averaging.  This  leads  to  the  'infinite'  resolution  of 
contact  discontinuities  and  shocks  that  the  scheme  displays. 

From  the  foregoing  discussion,  it  is  evident  that  the  success  of  the 
scheme  hinges,  to  a  large  extent,  on  the  inexpensive  and  exact  solution  of 
Rieraann  problems  and  an  appropriate  sampling  technique.   Ref.  (3)  describes  a 
modification  to  an  iterative  method  due  to  Godunov  (Ref.  5).   Theoretical 
details  for  the  Riemann  problem  solution  are  also  given  in  Ref.  (6). 

The  mathematical  properties  required  in  a  sampling  procedure  applicable 
to  this  scheme  are  defined  in  Ref.  (1).  A  brief  description  of  the  procedure 
is  given  below. 

In  previous  computations  using  the  RCM,  random  sampling  with  some 
variance  reduction  technique  (stratified  sampling);  was  used,  i.e.,  the  values 
were  taken  from  the  random  number  generator  installed  in  the  computer  (Ref. 
3).   It  was  shown  in  Ref.  (1)  that  a  more  accurate  form  of  sampling  is  a 
technique  due  to  van  der  Corput  (Ref.  7).   The  sequence  generated  is,  strictly 
speaking,  non-random,  but  has  particular  statistical  properties  that  are 
suitable  to  the  application.   The  sequence  is  referred  to  as  quasirandom  and 
is  generated  as  follows: 

The  binary  expansion  of  natural  numbers  may  be  expressed  as  (with   R=2): 
n  =  A0R°  +  A^1  +  A2R2  + +  AmR<n   ,   (0  v  A  k  <  R) 

ra 

I 
k=0 


i.e.    n  =  i     Ak.2k   ,  with  Ak  =  0  or  1   ,   n  =  1,  2,  3,  ... 


Next,    the   digits    of    the    binary   numbers    are    reversed   and   a   decimal   point    is    put 

preceding    the    number;    this    gives    the   numbers 

<Pn   =   AqR"1    +  AIR-2   +    ...    +  AnjR"^"1"1) 

ra 
or »  'hi   =  1     &k.2~(^+^)    ,  again  with  A^  =  0  or  1 
k=0 

Conversion  to  the  decimal  scale  of  these  numbers  yields  the  required  sequence 

of  quasirandom  numbers  defined  over  the  interval  [0,1],  i.e., 

<pn    (decimal)  =  0n  +  y 

or     0n  =  yn    (decimal)  -  y 

and    E,n   =  On.Ax  as  defined  earlier. 
The  first  few  elements  of  the  sequence  given  below  illustrate  the  construction 

n-1  (decimal)   =   1  (binary);   44   =   0.1  (binary)   =   U.5  (decimal) 

2  '  10  0.01  0.25 

3  11  0.11  0.7  5 

4  100  0.001  0.125 

5  101  0.101  0.625 

6  110  0.011  0.37  5 

7  111  0.111  0.875 

8  1000  0.0001  0.0625 


The  van  der  Corput  sequence  is  'equidistributed ' ,  and  yields  better  results 
than  those  obtained  using  a  'stratified'  random  sampling  technique. 

The  subroutine  employed  in  the  program  to  compute  the  random  numbers  is 
described  in  Appendix  B. 

2.2.   Boundary  Conditions 

In  general,  the  implementation  of  boundary  conditions  in  the  RCM  is 
quite  straightforward,  but  does  require  some  thought.   Referring  to  Fig.  2, 
the  b.c.'s  are  specified  at  points   1   and   N   for  the  left  and  right  boundary 


respectively.   Note  that  if  the  sampled  solution  at   (n+l)at   corresponds  to  a 
random  number   ^  <  0   ,  the  solution  is  assigned  to  the  grid  point  on  the 
left.   For  the  Riemann  problem  defined  by  points  1  and  3,  the  sampled  solution 
would  then  be  assigned  to  grid  point  1  at   (n+l)At   ;  however  this  is 
overridden  by  assigning  the  proper  boundary  condition  at  1  again,  and  there  is 
no  contradiction.   A  similar  procedure  is  adopted  at  the  right  hand  boundary 
when   Cn  >  0 

The  subroutines  for  the  boundary  conditions  are  named  in  the  format 
BCxn   ,   BC   standing  for  Boundary  Condition,   x  being  either  L  (for  Left), 
or  R  (for  Right  boundary)  and   n   being  a  number  from  I  to  5  with  the 
following  designations: 

1  -  solid  wall  condition 

2  -   outflow  at  constant  static  pressure 

3  -  special  formulation  ('piston'  inflow) 

4  -  isentropic  inflow  from  reservoir 

5  -  special  formulation  (rarefaction  wave  cancellation) 

2.2.1.   Solid  Wall  Conditions 

The  solid  wall  boundary  condition  requires  a  zero  normal  velocity  at 

the  wall  for  inviscid  flow  computations.   Due  to  the  random  sampling  involved 

Ax 
in  the  method  and  the  lateral  movement  of  the  sampled  solution  — —  to  the 

left  or  right  of  the  discontinuity,  the  condition  is  difficult  to  implement 
uniquely.   However,  the  procedure  adopted  here  is  found  to  yield  reasonably 
accurate  results  for  the  applications  intended.   (Note  that  the  difficulty  is 
not  unique  to  this  method  only.   The  implementation  of  zero  mass  flux  through 
a  surface  is  difficult  per  se  for  the  Euler  equations). 

Referring  to  Fig.  2,  let  the  physical  boundaries  be  at  point  2   and 


point  (N-l)  for  the  left  and  right  sides  respectively.   However,  the  boundary 
conditions  are  specified  at  point  1  (point  N)  for  the  left  (right)  side  as  a 
fictitious  'mirror'  state  of  the  conditions  at  point  3  (point  (n-2)) 
respectively,  but  with  the  reverse  sign  taken  for  the  velocity  component. 
Thus,  for  the  left  hand  boundary  Riemann  problem, 

PL  =  P(3)  ,  PL  =  P(3)  ,  uL  =  -u(3) 

PR  ■  P(3)  ,  PR  =  p(3)  ,  ur  =  u(3) 
and,  analogously,  for  the  right  hand  boundary  Riemann  problem. 

pL  =  p(N-2)  ,  mL  =  p(N-2)  ,  uL  =  u(N-2) 
pR  =  p(N-2)  ,  mr  =  p(N-2)  ,  uR  =  -u(N-2) 
The  solutions  are  then  sampled  in  the  manner  outlined  earlier. 

2.2.2.  Outflow  Conditions 

For  subsonic  outflow,  only  the  static  pressure   p   is  defined,  with 
the  continuation  condition  being  applied  to  the  rest  of  the  variables.   Thus, 
for  the  right  hand  boundary  for  example,  the  Riemann  problem  is  defined  as 
follows : 

PL  =  P(N~2)  ,  ml  =  P(N-2)  ,  uL  =  u(N-2) 
PR  =  Pout  >  PR  "  P(N"2)  ,  uR  =  -u(N-2) 
where   Pout   i-s  c^e  specified  outlet  pressure.   If  the  flow  going  out  is 
supersonic,  there  can  be  no  propagation  of  disturbances  upstream,  and  the 
continuation  condition  is  implemented  for  all  the  variables,  i.e.,  the  Riemann 
problem  now  is  the  trivial  case  defined  by 

PL  =  p(N-2)  ,  pL  -  p(N-2)  ,  uL  =  u(N-2) 
Pr  =  p(N-2)  ,  pr  =  p(N-2)  ,  uR  =  -u(N-2) 

2.2.3.  Special  Formulation  of  'Piston'  Inflow 

In  general,  for  idealized  wave  rotor  flows,  hot  combustion  gases  are 
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introduced  into  the  rotor  through  nozzles  angled  such  as  to  allow  the  flow  to 
'slip  onto'  the  rotor,  i.e.,  without  incurring  incidence  or  deviation  angle 
losses.   Also,  in  the  ideal  treatment,  the  air  in  the  passages  of  a  wave  rotor 
is  exposed  to  the  hot  gas  at  high  pressure  instantaneously.   The  idealizations 
allow  for  uniform  conditions  to  be  prescribed  at  the  hot  gas  inlet  port. 
Thus,  a  'special'  form  of  inflow  boundary  condition  needs  to  be  specified 
here,  namely,  the  static  pressure,  the  velocity  and  the  density  of  the 
incoming  hot  gas.   Although  equivalent  to  specifying  the  total  pressure  and 
temperature  in  the  usual  inflow  boundary  condition  treatment,  some  thought  is 
required  in  wave  rotor  type  flows  when  specifying  pJas   ,   pgas   an^   naas 
This  is  because  only  a  shock  wave  needs  to  be  generated,  with  no  waves 
travelling  opposite  to  the  direction  of  flow.   The  solution  to  the  Rieniann 
problem  would  then  consist  of  just  two  states  connected  by  a  single  shock 
wave.   The  flow  is  equivalent  to  that  generated  when  a  piston  is  pushed 
instantaneously  into  a  gas  at  rest.   In  general,  the  state  of  the  air  inside 
the  rotor  passage  is  known;  explicit  relations  for  two  states  connected 
through  a  shock  wave  are  given  in  Ref.  (6).   These  so-called  transition 
functions  help  in  specifying  the  boundary  conditions  for  the  incoming  flow 
properly. 

If  we  consider  the  left  boundary  for  this  inflow,  the  Rieraann  problem 
is  set  up  as: 

PL  =  Phot  gas  »  ^L  =  Phot  gas  >  UL  =  uhot  gas 
PR  ■  P(3)  ,  PR  -  p(3)  ,  ur  =  u(3) 
with   pl   ,   pl  and  u^   having  been  chosen  in  accordance  with  the 
considerations  discussed  above. 
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2.2.4.   Isentropic  Inflow  From  Reservoir 

The  induction  of  fresh  charge  or  air  onto  the  rotor  usually 
corresponds  to  an  isentropic  inflow  situation.   The  flow  in  the  vicinity  of 
the  passage  end  can  be  treated  as  quasi-steady,  with  the  assumption  that  no 
flow  separation  takes  place  when  the  flow  enters.   Two  boundary  conditions  are 
required  for  this  type  of  inflow;  these  are  provided  by  the  conservation  of 
energy  in  the  flow  from  the  external  region  to  the  inlet  (assumed  to  be 
steady),  and  by  the  prescibed  entropy  level  of  the  gas  in  the  external 
region. 

The  boundary  conditions  may  thus  be  expressed  as 

2     2    2       2    2 
u.   +  — p  a,    =  — r  a 
in    y-1   in     Y~l   tot 

sin  =  stot 
where  the  subscripts  'in'  and  'tot1  apply  to  conditions  at  the  inlet  of  the 

passage  and  external  reservoir  respectively.   The  sonic  velocity  is  denoted  by 

a   ,  and  flow  velocity  by  u   .   Note  that  knowledge  of  the  Riemann  variable 

arriving  at  the  passage  end  from  within  the  passage  is  required  to  be  able  to 

solve  the  energy  equation  above  for   a^n   and  uin   .   For  the  left  end,  for 

example , 

2 
Qin  "  -yTf  ain  ~  uin 

which  together  with  the  energy  equations  yields 


ai 


J  Y+l   2     y-1   2 
n  -  Qin  +  ^  Y-i  atot     2   QLn 


Y+l 
Y-1 

and  subsequently  the  other  variables. 
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The  simple  analytical  treatment  given  above  has  to  be  modified 
somewhat  if  a  contact  discontinuity  is  formed  when  the  inflow  begins.   This 
is  due  to  the  fact  that  the  value  of  the  arriving  Riemann  variable  is  changed 
across  such  a  discontinuity,  which  thus  leads  to  an  additional  unknown. 
Procedures  for  solving  the  inflow  for  these  situations  are  given  in  Ref.  (8). 
In  the  program  developed  here,  reasonably  good  results  are  obtained  by  setting 
the  velocity  at  the  boundary  point  equal  to  the  velocity  at  the  point  nearest 
the  physical  boundary.   For  the  left  end  e.g.,  the  variables  for  the  Left 
state  of  the  Riemann  problem  are  obtained  as  follows: 

u(l)  =  u(3)   , 
a  reasonably  accurate  assumption  just  at  the  point  of  inlet  opening. 

Then,  from  the  'energy  ellipse', 


a<1)  "  Vatot   -  *£     u(l)2 


MO)  =    )  (    »  incoming  Mach  number 
a  (  1 ) 

n(\)    =  Ptot 

P  Y-l       9   Y/(Y-1)  ' 

[1+Y-  M(l)2]  y/u  i; 

with  similar  isentropic  relations  to  compute  other  flow  variables.   Note 
that  once  the  interface  or  contact  discontinuity  has  moved  a  certain  distance 
inside  the  passage,  the  simple  analytical  expressions  given  earlier  in  the 
section  can  be  used,  since  now  the  value  of  the  arriving  Riemann  variable 
would  be  known  at  the  boundary. 
2.2.5.   Special  Formulation  for  Rarefaction  Wave  Cancellation 

The  spreading  of  rarefaction  fans  leads  to  unwanted  wave  reflections 
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which  occupy  large  zones  in  the  passages  of  wave  rotors.   Fig.  (4)  shows  a 
wave  diagram  proposed  by  Spectra  Technology,  Inc.,  which  incorporates 
so-called  'wave  management'  or  'tuning'  ports  to  ideally  cancel  (and  otherwise 
attenuate)  impinging  rarefaction  fans.   The  physical  boundary  conditions  are 
thus  dictated  by  the  flow  developing  in  the  passage,  i.e.,  the  port  has 
non-uniform  flow  conditions  in  it,  which  at  each  point  match  those  of  the  flow 
at  the  end  of  the  passage  so  as  to  disallow  any  reflections  to  take  place. 
Numerically,  this  condition  is  achieved  by  implementing  the  continuity 
condition  across  the  boundary  for  all  the  flow  variables  involved.   For  the 
left  boundary,  thus,  the  Riemann  problem  is  defined  by: 

PL  =  p(3)   ,   PL  =  P(3)   ,   ul  =  u(3) 

PR  =  p(3)   ,   PR  =  P(3)   ,   uR  =  u(3) 
and  analogously  for  the  right  boundary.   Note  that  these  boundary  conditions 
may  involve  either  inflow  or  outflow. 

2.3   Example  Calculations 

The  listing  of  the  program  is  included  in  Appendix  A,  and  the  various 
names  for  the  variables  are  listed  in  Appendix  B,  along  with  some  instructions 
on  how  to  use  the  program.   No  effort  as  yet  has  been  made  to  optimize  the 
code  either  for  storage  requirements  or  for  execution  efficiency. 

In  this  section,  some  sample  calculations  are  carried  out  using  the  code, 
to  illustrate  its  usefulness  in  constructing  idealized  design  point  wave 
diagrams  which  can  serve  as  the  starting  configuration  for  detailed 
construction  of  diagrams  incorporating  real  flow  effects. 

2.3.1.   Test  Case  for  1-D,  Inviscld,  Unsteady,  Compressible  Flow 

Fig.  (1)  illustrates  the  initial  conditions  in  a  shock  tube,  with  the 
diaphragm  at   xq   .   Sod  (Ref.  9)  suggested  a  test  case  for  hyperbolic 
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conservation  laws  with  the  following  conditions  as  initial  states  in  the  shock 
tube : 

PI  =  1.0   ,   pi  =  1.0   ,  ui  =  0.0 
P5  =  0.1   ,   P5  =  0.125   ,   U5  =  0.0 
i.e.,  the  gas  on  either  side  of  the  diaphragm  is  in  a  quiescent  state 
initially.   The  ratio  of  specific  heats  is  chosen  to  be  7/5,  and   Ax   is 
chosen  to  be  0.01. 

The  solution  (before  any  wave  has  reached  either  the  left  or  right 
end)  is  shown  in  Fig.  (5).   The  squares  shown  at  locations   x^   ,   x2   ,   X3 
and   x/,   in  the  density  plot  give  the  analytically  calculated  amplitude  and 
location  of  the  head  -  and  tail  waves  of  the  left-running  r.iref  action ,  the 
contact  surface  moving  to  the  right  and  the  shock,  wave  moving  at  supersonic 
velocity  to  the  right  respectively.   The  solid  lines  are  the  solutions 
obtained  by  the  RCM  at  different  time  levels;  the  zero  numerical  diffusion 
feature  of  the  method  is  evident  in  the  'infinite'  resolution  of  the  contact 
discontinuity  and  the  shock,  and  the  dispersion  (phase  error)  is  within  one 
grid  spacing.   The  constant  states  are  perfectly  realized. 

It  is  these  features  of  the  method  that  make  it  very  attractive  for 
application  to  wave  rotor  type  flows,  since  the  successful  design  of  the 
device  is  predicated  on  being  able  to  accurately  compute  wave  arrival  times  at 
the  various  ports. 

2.3.2.   Wave  Turbine  Experiment 

Ref.  (10)  describes  the  wave  rotor  experimental  set  up  at  the 
Turbopropulsion  Laboratory.   Initial  tests  being  carried  out  currently  are 
with  the  wave  rotor  in  a  turbine  mode,  i.e.,  one  side  of  the  rotor  is  blocked 
off,  and  high  pressure  air  is  brought  onto  the  rotor  and  taken  off  again  from 
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the  other  side.   The  passages  of  the  rotor  being  angled  at  60°  to  the  axis, 
the  180°  reversal  in  the  direction  of  the  fluid  flow  creates  an  angular 
momentum  change,  in  turn  generating  large  turbomachinery  work,  coefficients. 
Fig.  (6)  shows  the  wave  diagram  computed  using  the  code.   The  movement  of  the 
rotor  is  from  top  to  bottom.   At   t=0   ,  the  high  pressure  air  is  brought 
into  contact  with  quiescent  atmospheric  air  in  the  rotor  passages,  at  point  a. 
This  corresponds  to  the  'piston'  inflow  boundary  condition  described  in 
section  2.2.3..   A  shock,   S   ,  is  generated  immediately,  (idealized  case  of 
instantaneous  cell  opening),  which  travels  from  the  right  to  the  left,  and 
strikes  the  solid  wall  at  the  left  end.   The  reflection  of  the  shock  takes 
place  at  point  b  according  to  the  solid  wall  boundary  condition  described  in 
section  2.2.1..   Behind  this  shock,  and  moving  at  a  slower  velocity  is  the 
contact  surface,   T   ,  which  penetrates  into  the  passage  only  a  fractional 
distance  before  encountering  the  reflected  shock,  RS ,  at  point  c.   The 
reflected  shock  is  transmitted  through  the  contact  surface,  (bringing  the  flow 
to  a  near  halt),  and  reaches  the  right  side  at  point  d,  whereupon  the  inlet 
port  is  closed.   The  air  trapped  in  the  rotor  passages  is  now  at  a  high 
pressure  and  in  a  quiescent  state.   When  this  air  is  released  at  point  e  to  a 
low  pressure  region,  a  rarefaction  wave  is  generated,  R  ,  which  travels  to  the 
left,  spreading  out  in  the  process.   It  interacts  with  the  stationary  contact 
surface,  I  ,  setting  it  into  motion  again,  and  reflects  off  the  solid  wall  at 
the  left  as   RR   .   The  boundary  condition  imposed  at  point   e   is  the  outflow 
at  constant  static  pressure  condition  described  in  section  2.2.2..   The  outlet 
port  is  closed  at  a  time  when  the  exit  velocity  falls  to  about  half  its 
initial  value. 

This  experiment  embodies  two  fundamental  processes  in  wave  rotors: 
those  of  cell  filling  and  cell  emptying.   Almost  all  the  other  processes 
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typical  to  wave  rotors  are  combinations  of  the  cell  filling  and  cell  emptying 
unit  processes.   Comparison  of  the  ideal  computed  numbers  obtained  here  with 
experimental  data  will  provide  information  helpful  in  the  identification  and 
sources  of  losses. 

The  program  is  set  up  to  start  at   t=0   in  this  case,  with  initial 
data  provided  along  the  entire  passage,  i.e.,  from  x=0   to   x=0. 1863m   (the 
actual  length  of  the  wave  rotor  being  tested).   Since  the  passages  have 
quiescent  atmospheric  air  in  them  at   t=0   ,  the  initial  data,  of  course, 
describes  these  conditions.   Switches  for  the  left  and  right  boundaries 
describe  what  type  of  boundary  conditions  prevail  and  direct  the  program  to 
the  appropriate  subroutines.   These  switches,  designated   SWL   and   SWR   ,  for 
left  and  right  respectively,  are  assigned  integer  number  values  which 
correspond  to  the  numeric  value  of  the  particular  boundary  condition  they 
represent.   Thus,  if  the  left  boundary  is  a  solid  wall,   SWL= 1   , 
corresponding  to  the  boundary  condition  subroutine   BCL1   .   In  this  example 
then,  the  initial  switch  settings  at   t=0   are   SWL=1   and   SWR=3   , 
corresponding  to  a  solid  wall  at  the  left  and  a  'piston'  inflow  at  the  right 
(which  starts  at   t=0  at  point   a).   At  point   d   ,  the  switches  are  reset  to 
SWL-1   and   SWR=1   due  to  the  closure  of  the  inlet  port.   At  point   e   ,  the 
switches  are   SWL=1   ,   SWR=2   ,  signifying  opening  of  the  exhaust  port  with 
outflow  at  a  constant  static  pressure.   The  whole  wave  diagram  can  thus  be 
packaged  into  a  'module'  subroutine  and  called  from  the  main  program  with  a 
single  call  statement.   This  type  of  modularity  allows  for  wave  diagrams  of 
different  'families'  to  be  developed  by  simply  calling  the  right  'module' 
subroutine. 

The  next  two  examples  illustrate  this  concept  as  they  deal  with  two 
very  different  types  of  wave  diagrams. 
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2.3.3.   General  Electric  Wave  Engine 

Fig.  (7)  shows  a  schematic  of  the  wave  diagram  constructed  for  the 
G.E.  wave  engine.   Briefly,  the  device  is  configured  for  a  gas  generator  mode 
of  operation,  with  counterflow  scavenging,  and  is  capable  of  producing  net 
shaft  power.   For  a  fuller  description  of  the  machine,  see  Ref.  (11).   In  this 
example,  fresh  charge  (air)  is  induced  into  the  rotor  (from  an  external 
reservoir)  through  the  wave  action  of  the  rarefaction  fan  originating  at  the 
exhaust  port  opening.   The  usefulness  of  the  rotor  is  gauged  by  the  net 
pressure  rise  across  the  machine,  i.e.,  the  ratio  of  the  total  exhaust 
pressure  to  total  (fresh  air)  inlet  pressure. 

For  performance  estimation  purposes,  it  is  sufficient  to  investigate 
only  the  exhaust  and  induction  processes  as  shown  in  Fig.  (8).   The  initial 
data  specified  is  as  follows:   the  exhausting  pressure  ratio   Pe/Po   >  tne 
total  pressure  ratio  across  the  rotor   Pte/Pta   and  an  assumed  total 
temperature  ratio  Tte/Tta   •   *n  tni-s  particular  cycle,  the  amount  of  fresh 
charge  induced  in  is  ideally  equal  to  the  gases  exhausted  out,  i.e., 
min  =  "but   >  and  this  mass  balance  is  carried  out  after  each  computation  to 
correct  the  assumed  temperature  ratio  Tte/Tta   (which  otherwise  constitutes 
overspecif ication  of  the  initial  conditions). 

The  calculation  starts  at   t=0   ,  with  initial  data  consistent  with 
the  chosen  pressure  and  temperature  ratios  specified  along  the  passage  length. 
Initial  switch  settings  are   SWL=1   and   SWR=2   for  the  solid  wall  boundary  at 
the  left  and  the  exhaust  to  a  constant  pressure  at  the  right.   As  shown  in  the 
figure,  a  rarefaction  fan  is  generated,  propagating  to  the  left  and  reflecting 
off  the  solid  wall.   At  time   t=Tj   >  tne  pressure  at  the  wall  has  been 
reduced  to  that  outside  the  passage,   pta   ,  which  is  when  the  inlet  port  is 
opened.   The  switches  are  now  set  to   SWL=4   and   SWR=2   for  isentropic  inflow 
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from  an  external  reservoir  at  the  left,  and  still  outflow  at  a  constant 
pressure  at  the  right.   The  exhaust  port  is  closed  at  time   t=T2  which 
corresponds  to  the  exit  velocity  having  dropped  off  to  approximately  half  its 
steady  state  value  at  the  beginning  of  the  exhaust  process.   Now  the  switches 
are  set  to   SWL=4   and   SWR=1   ,  for  the  solid  wall  condition  at  the  right. 
The  sudden  closure  of  the  exhaust  port  generates  a  'hammer'  shock  travelling 
to  the  left,  interacting  with  the  Incoming  interface  (shown  by  dashed  line), 
and  reaching  the  passage  end  at   t=T4   at  which  time  the  inlet  port  is  closed, 
with  the  switches  being  reset  to   SWL= 1   and   SWR=1   .   Note  the  reflected 
shock  travelling  from  left  to  right  generated  at  the  interaction  of  the 
contact  surface  and  the  hammer  shock. 

Once  this  solution  is  obtained,  integration  of  the  mass  flux  through 
the  inlet  and  exhaust  ports  is  carried  out  and  if  the  two  numbers  do  not 
match,  the  assumed  temperature  ratio  Ite/^ta   *s  adjusted  in  the  initial 
data,  till  such  time  as  m±n   =   mOK1^ 

This  calculation  is  sufficient  for  performance  analyses:   If  the 
entire  wave  diagram  has  to  be  worked  out,  then  at  a  time   t  >  C3   ,  hot  gas 
from  the  combustion  chamber  is  brought  onto  the  rotor  (the  boundary  condition 
corresponding  to  'piston'  inflow)  on  the  right  hand  side.   This  would 
generate  the  shock  to  compress  the  induced  air  and  when  this  shock  reached 
the  left  end,  the  transfer  port  (see  Fig.  7)  would  be  opened  for  such  time  it 
takes  for  the  compressed  air  to  be  completely  scavenged  out  of  the  rotor. 
Fig.  (9)  shows  some  performance  curves  obtained  using  the  procedure  outlined 
above.   In  Figs.  (10a,  b,  c)  are  shown  three  sets  of  flow  parameters  at 
different  time  steps  corresponding  to  the  inlet  port  just  opening,  the 
exhaust  port  closing  and  the  inlet  port  closing;  the  qualitative  distributions 
of  the  flow  parameters  in  the  passage  are  immediately  seen  to  be  accurate  when 
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compared  with  the  wave  diagram  shown  in  Fig.  (8).   Of  interest  is  the  set  of 
plots  for  the  time  step  when  the  inlet  port  has  just  been  closed.   The  flow 
between  the  end  of  the  passage  and  the  location  of  the  interface  is  seen  to 
be  quite  non-uniform  in  the  density  plot.   At  the  same  time,  the  shock 
reflected  from  the  interface  has  reached  the  right  side  and  reflected  off  the 
solid  wall.   These  considerations  help  to  decide  optimum  port  opening  and 
closing  times.   For  example,  Fig.  (11)  shows  what  happens  if  the  inlet  port  is 
not  closed  at  just  the  time  the  shock  reaches  the  end,  but  rather  at  some 
short  time  later.   The  shock  now  sees  an  open  boundary  and  reflects  off  as  an 
expansion  to  match  the  high  pressure  behind  it  with  the  incoming  total 
pressure  which  is  at  a  lower  value.   This  reflected  expansion  is  manifested  in 
the  pressure,  density  and  velocity  plots  of  the  figure. 

The  Entire  sequence  of  wave  interactions  of  this  example  is  computed 
by  the  RCM  without  the  implementation  of  artificial  viscosity  or  artificial 
compression  methods,  or  tracking  and  capturing  schemes.   This  'hands  off' 
feature  of  the  method  renders  it  eminently  useful  for  fast  preliminary 
evaluations  of  complex  wave  diagrams  for  the  application  at  hand. 

The  next  example  computes  an  idealized  wave  diagram  for  the  nine-port 
pressure  exchanger  concept  proposed  by  Spectra  Technology,  Ref.  (12). 

2.3.4.   Spectra  Technology  Pressure  Exchanger 

Fig.  (4)  shows  the  ideal  wave  diagram  for  the  nine-port  pressure 
exchanger.   This  configuration  is  a  good  case  example  to  compute  with  the  RCM 
because  of  the  different  types  of  boundary  conditions  that  need  to  be  dealt 
with  in  the  evaluation  of  the  cycle.   The  computation  is  started  at   t=0   ,  at 
the  point  of  high  pressure  hot  gas  inlet  (driver  gas  inlet).   In  the  manner 
described  in  the  G.E.  wave  engine  example,  the  initial  data  is  prescribed  for 
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the  entire  passage  at  this  time  step  and  the  boundary  condition  switches  are 
initially  set  at   SWL=1   and   SWR=3   for  the  solid  wall  at  the  left,  and  the 
'piston'  inflow  at  the  right  hand  end.   Since  there  is  a  multiplicity  of  types 
of  boundary  conditions,  e.g.,  three  outflow  ports,  an  index,  JCOUNT,  is  used 
to  ensure  proper  sequencing  of  the  switches.   The  following  table  is 
presented  as  an  example  of  the  settings  of  the  switches  to  carry  out 
calculations  for  one  cycle.   The  inflow  and  outflow  port  conditions  are  those 
proposed  by  Spectra  Technology  for  their  idealized  diagram. 


TIME  STEP,  N 

JCOUNT 

SWL 

SWR 

REMARKS 

0 

0 

1 

3 

CYCLE  STARTS.   HP  GAS  INLET  PORT  OPENS 

500 

1 

2 

3 

HP  AIR  OUTLET  PORT  OPENS 

1408 

2 

2 

1 

HP  GAS  INLET  PORT  CLOSES 

1765 

3 

5 

1 

HP  AIR  OUTLET  PORT  CLOSES.   TUNING  PORT 
LI  OPENS 

1816 

4 

2 

1 

TUNING  PORT  LI  CLOSES.   IP  GAS  OUTLET 
PORT  OPENS  (PORT  El) 

2069 

5 

2 

5 

TUNING  PORT  Rl  OPENS 

2261 

6 

2 

I 

TUNING  PORT  Rl  CLOSES 

2595 

7 

5 

1 

IP  GAS  OUTLET  PORT  CLOSES.   TUNING  PORT 
L2  OPENS 

2636 

8 

2 

1 

TUNING  PORT  L2  CLOSES.   LP  GAS  OUTLET 
PORT  OPENS  (PORT  E2) 

3029 

9 

2 

5 

TUNING  PORT  K2  OPENS 

3237 

10 

2 

4 

TUNING  PORT  R2  CLOSES.   LP  AIR  INLET 
PORT  OPENS 

4961 

11 

1 

4 

LP  GAS  OUTLET  PORT  CLOSES 

5529,0 

0 

1 

3 

LP  AIR  INLET  PORT  CLOSES.   CYCLE 
COMPLETED 

The  total  cycle  time  as  calculated  by  the  RCM  is  3.0676  mseconds,  which 
compares  well  with  the  time  computed  by  Spectra  Technology  (using  the 
FCT-SHASTA  algorithm)  of  3.07  mseconds.   The  execution  time  on  an  IBM 
370-3033AP  for  the  5529  steps  computed  in  the  example  above  was  3  minutes 
38  seconds,  including  the  I/O  operations  and  the  graphics. 


21 


Figs.  (12a,  b,  c)  show  three  sets  of  plots  of  the  flow  parameters  for 
the  following  cases:   a)  the  H.P.  air  outlet  port  opens  on  time,  i.e.,  just  as 
the  shock  reaches  the  left  end  of  the  passage,  b)  the  port  opens  before  the 
shock  has  reached  the  end,  and  c)  the  port  opens  after  the  shock  has  reached 
the  end.   The  constant  pressure  and  velocity  states  that  prevail  in  the 
passage  just  after  the  shock  has  reached  the  left  end  (time  'section'  line   tj 
in  wave  diagram),  are  perfectly  realized  in  Fig.  (12a),  while  the  contact 
surface  is  at  the  location  shown  by  the  sharp  discontinuities  in  the  density 
and  entropy  plots.   Should  the  inlet  port  be  opened  earlier,  e.g.,  at  the  time 
level  shown  by   x^_   in  the  wave  diagram,  what  happens  is  as  follows:   the 
pressure  in  the  passage  is  still  at  the  pre-corapressed  level  and  this  comes 
into  contact  with  the  pressure  level  in  the  port  which  is  considerably  higher, 
resulting  in' a  shock  propagating  into  the  passage,  colliding  with  the  left 
moving  shock  and  raising  the  overall  pressure  level  to   ~3.0  as  shown  in 
Fig.  (12b).   However,  as  soon  as  the  left  moving  shock  readies  the  end,  it  now 
encounters  an  open  boundary  with  conditions  that  do  not  match  those  behind  the 
shock,  resulting  in  a  rarefaction  fan  being  generated,  which  propagates  to  the 
right.   This  expansion  fan,  travelling  at  sonic  velocity  relative  to  the  gas 
into  which  it  is  propagating,  soon  overtakes  the  right  moving  shock  which  is 
travelling  at  a  subsonic  velocity  relative  to  the  same  gas.   This  interaction 
results  in  an  attenuation  of  both  the  rarefaction  as  well  as  the  shock  wave. 
Note  that  the  overall  pressure  and  velocity  levels  behind  the  rarefaction  are 
about  the  same  as  for  case  a),  i.e.,  the  effects  of  the  mismatch  are  not  very 
significant  at  the  outlet  port.   However,  should  the  right  moving  pressure 
perturbations  of  case  b)  not  attenuate  each  other  significantly  before  they 
reach  the  right  hand  end,  the  consequences  could  be  severe  for  the  overall 
wave  diagram,  since  this  will  lead  to  further  (unwanted)  wave  reflections. 
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Fig.  (12c)  shows  what  occurs  if  the  outlet  port  is  opened  too  late, 
corresponding  to  time  level  t\+     on  the  wave  diagram.   Now  the  left 
travelling  shock,  encounters  a  wall  boundary  condition  on  reaching  the  left 
end  and  reflects  off  as  a  shock,  effectively  doubling  the  pressure  level 
behind  it  (  >3.5  in  pressure  plot  of  Fig.  (12c)).   When  the  outlet  port  opens, 
there  is  again  a  mismatch  of  conditions  in  the  port  and  in  the  passage,  with 
the  pressure  level  in  the  passage  being  considerably  higher  than  that 
prescribed  for  the  outlet  port.   A  rarefaction  wave  is  generated  which 
propagates  to  the  right  and  overtakes  the  reflected  shock.   The  same  criterion 
holds  for  this  case  too,  i.e.,  the  ensuing  attenuation  of  these  pressure 
pulses  should  occur  before  they  reach  the  right  hand  end,  preferably  even 
before  they  reach  the  interface  still  propagating  towards  the  left  at  the 
flow  velocity. 

The  considerations  above  give  a  preview  of  the  nature  of  decisions 
required  in  the  successful  design  of  a  wave  rotor  device.   It  is  clear  that 
quite  a  few  iterations  are  involved  in  the  process  of  designing  a  viable  wave 
diagram  for  a  particular  application,  and  each  iteration  entails  calculating 
two  or  more  complete  cycles  to  ensure  'closure'  or  repeatability  of  the  cycle. 
A  fast  solver  like  the  RCM  allows  reaching  an  idealized  'base'  design  quickly 
and  inexpensively. 

Appendix  A  is  a  listing  of  the  program  in  its  present  development 
stage.   As  mentioned  earlier,  no  attempt  has  been  made  to  optimize  the 
program,  either  for  storage  requirements  or  for  execution. 

Appendix  B  gives  a  description  of  the  structure  of  the  program,  a 

listing  of  the  important  variables,  the  subroutines  and  the  function 

subprograms.   A  step  by  step  guide  is  also  included  to  set  up  and  run  the 

program. 
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3.   DISCUSSION  AND  RECOMMENDATIONS 

3.1.  Discussion 

For  meeting  the  criteria  listed  in  the  Introduction,  in  one  dimension, 
Glimm's  method  or  the  RCM  appears  to  be  superior  to  any  difference  method. 
For  wave  rotor  type  applications,  where  discontinuities  need  to  be  computed 
with  sharpness,  the  'infinite'  resolution  of  such  discontinuities  inherent  in 
the  RCM  make  it  a  natural  choice  to  carry  out  ideal  flow  calculations  for 
preliminary  design  purposes.   Boundary  conditions  can  be  implemented  quite 
easily  and  do  not  require  information  from  points  outside  the  domain  of 
dependance  as  is  the  case  in  some  finite  difference  schemes.   The  van  der 
Corput  sampling  technique  results  in  the  best  possible  representation  of  the 
wave  propagation,  which  is  essential  for  the  correct  representation  of 
continuous  waves,  particularly  those  produced  by  nonlinear  interactions. 

The  method,  however,  is  not  recommended  to  solve  for  f Lows  with  real 
effects  such  as  friction,  heat  transfer  and  area  change,  or  to  be  extended  to 
multi-di mensional  flows.   Although  considerable  research  is  being  done  to 
rigorously  extend  the  method  to  such  flows,  with  some  degree  of  success  (see 
Refs.  1,  4,  13),  the  present  state  of  development  is  not  mature  enough  to 
ensure  a  useful  practical  code  as  the  outcome. 

3.2.  Recommendations 

Many  options  are  available  for  one  wishing  to  develop  either  a  1-D  code 
with  real  effects  and/or  a  multi-dimensional  code  for  wave  rotor  type 
applications.   The  author  prefers  to  recommend  numerical  formulations  which 
are  dependent  on  the  solution  of  Riemann  problems,  such  as  the  Godunov  method; 
the  motivating  reason  for  this  preference  is  that  a  Riemann  problem 
constitutes  the  solution  of  a  discontinuity  in  the  flow  in  terms  of  other 
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discontinuities  (if  any  are,  indeed,  present),  and  the  scheme  is  thus 
intrinsically  suited  for  solving  such  flows;  on  the  other  hand,  the  other 
schemes,  in  general,  require  to  be  made  aware  of  discontinuities  in  the  flow 
through  some  external  device,  and  then  treat  them  through  other  artificial 
devices. 

A  second-order,  quasi  one-dimensional  (variable  cross-sectional  area) 
scheme  has  recently  been  developed  by  Ben-Artzi  and  Falcovitz  (Ref.  14).   The 
method  is  based  on  the  exact  solution  of  'generalized  Rieraann  propbleras ' ,  and 
has  demonstrated  very  good  results;  it's  least  accurate  approximation  is 
equivalent  to  Godunov's  first  order  method  (Ref.  9).   The  resolution  of  shocks 
and  other  discontinuities  and  singularities  of  the  flow  field  is  also  high. 
Extension  to  more  than  one  dimension  appears  to  be  straightforward  through  the 
use  of  operator  splitting  techniques,  but  has  as  yet  not  been  tried 
extensively. 
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Fig. 8  :  Gas  Exliaust  and  Fresh  Air  Induction 
Process  in  C.E.    Wave  Engine 
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Fig. 12  a  :  Distribution  of  Flow  Parameters  in 
Rotor  Passage  when  the  H.P.  Air 
Outlet  Port  Opens  on  Time. 
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APPENDIX  A 


Listing  of  Program  RCM 


40 


PROGRAM  RCM  WITH  VAN  DER  CORPUT  SAMPLING  AND  SINGLE  TIME  STEP  RCM00030 

INTEGER  QPRINT,  QSTOP ,  SWL,  SWR  RCM00040 

DIMENSION  XX(6) ,YY( 6)  RCM00050 

DIMENSION  XARRAY(IOO)  RCM00060 

DIMENSION  WNORM(12),IDIGT(12)  RCM00070 

DIMENSION  P(203),R(203),U(203),A(203),S(203),X(203)  RCM00080 

COMMON/SUBS/P,R,U,A,S,X  RCM00090 
COMMON/ GLIMM1/ PGLIM , RGLIM , UGLIM , PL , RL , UL , PR , RR , UR , AL , AR , GL , GR , EPS  RCMOO 100 

COMMON /GLIMM2/DT,DX, XI  RCM00110 

COMMON/ FUN1/ G , PA , RA , UA , RB , RMU  RCMOO 120 

COMMON/ S AMPLE /WNORM,IDIGT  RCM00130 

COMMON  XARRAY,N1  RCMOO 140 

CALL  COMPRS  RCMOO 150 

CALL  BLOWUP (0.5)  RCM00160 

CALL  PAGE (11. 0,8. 5)  RCMOO 170 

CALL  HWSCAL( 'SCREEN' )  RCM00180 

DATA  K,SWL,SWR/500,1,3/  RCM00190 

DATA  N,CFLNUM,TTOTAL/ 0,0. 60,0.0/  RCM00200 
DATA  PSEXIT,PSINL,PSINR,RINL,RINR/116954. ,3770000. , 3819952 . 50 , 6 . 8 ,RCM00210 

--6.800/  RCM00220 

DATA  PSOUTl,PSOUT2,PSOUT3/ 38 19952. 5, 243 1800. 0,15 30007. 5/  RCM00230 

DATA  PTOTIN,RTOTIN/1656663.8,7.486/  RCM00240 

DATA  PREF,RREF,XREF/1656663.8,7.486,0. 1800/  RCM00250 

G=1.4  RCM00260 

GL=1.4  RCM00270 

GR=1.4  RCM00280 

EPS=l.E-06  RCM0O290 

QSTOP=20  RCM00300 

N1  =  0  RCM00310 

JCOUNT=0  RCM00320 

KCOUNT=0  RCM00330 

UEXMAX=0.  RCM00340 

DX=0.01  RCM00350 

AREF=SQRT(PREF/RREF)  RCM00360 

TIMREF=XREF/AREF  RCM00370 

RMU=SQRT((G-1.)/(G+1.))  RCM00380 

X(l)=-0.5-DX  RCM00390 

ZETA=WDP(1)  RCM00400 

XII=DX*(WDP(0)-0.5)  RCM00410 

DO  25  1=2,203  RCM00420 

X(I)=X(I-1)+0.5*DX  RCM00430 

25  CONTINUE  RCM00440 

DO  35  1=1,100  RCM00450 

XARRAY(I)=X(I*2+1)  RCM00460 

35  CONTINUE  RCM00470 

INITIAL  DATA  RCM00480 

CALL  INIT1  RCM00490 

CALL  INIT2L(PSEXIT)  RCM00500 

CALL  INIT2R(PSEXIT,PREF,RREF)  RCM00510 

CALL  INIT3L(PSINL,RINL)  RCM00520 

CALL  INIT3R(PSINR,RINR)  RCM0O530 

NONDIMENSIONALIZATION  RCM00540 

DO  30  1=1,203,2  RCM00550 

P(I)=P(I)/PREF  RCM00560 
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R(I)=R(I)/RREF 

U(I)=U(I)/AREF 

A(I)=A(I)/AREF 

S ( I ) =  ALOG ( P ( I ) /R ( I ) * *G ) 
30  CONTINUE 

CALL  PLOTl(K) 

DO  40  J=1,K 

N  =  N+1 

XII=DX*(WDP(0)-0.5) 

QPRINT=N/50 

DT=100. 

DO  50  1=1,203,2 

DTT=CFLNUM-DX/(2."AMAX1(ABS(A(I)+U(I)),ABS(A(I)-U(I)))) 

DT=AMIN1(DTT,DT) 
50  CONTINUE 

TTOTAL=TTOTAL+DT 

TIME=TTOTAL*TIMREF 

XI=-XII 

DO  60  1=1,201,2 

PL=P(I) 

RL=R(I) 

UL=U(I) 

PR=P(I+2) 

RR=R(I+2) 

UR=U(I+2) 

XITEMP=XI 

IF(I.EQ.l)  XI=ABS(XI) 

IF((I.EQ.201).AND. (XI.GT.0.0))  XI=-XI 

CALL  GLIMM(QSTOP,PSTAR,USTAR,ASTAR) 

XI=XITEMP 

P(I+1)=PGLIM 

R(I+1)=RGLIM 

U(I+1)=UGLIM 
60  CONTINUE 

DO  70  1=1,201,2 

IF(XI.LT.O. )  GOTO  80 

P(I+2)=P(I+1) 

R(I+2)=R(I+1) 

U(I+2)=U(I+1) 

A(I+2)=SQRT(G*P(I+2)/R(I+2)) 

S(I+2)=ALOG(P(I+2)/R(I+2)**G) 

GOTO  70 
80  P(I)=P(I+1) 

R(I)=R(I+1) 

U(I)=U(I+1) 

A(I)=SQRT(G»P(I)/R(I)) 

S(I)=ALOG(P(I)/R(I)**G) 
70  CONTINUE 
C      CALL  GE(SWL,SWR,N,TTOTAL,TIME,UEXMAX,PTOTIN,PREF) 
C      CALL  DETON(SWL,SWR,N,QPRINT,TTOTAL,TIME) 

CALL  SPCTRA (N , SWL , SWR , TIME , UEXMAX , PSEXIT , PSOUT1 , PSOUT2 , PSOUT3 , J 

-COUNT , QPRINT , TTOTAL , KCOUNT ) 

IF(SWL.EQ.i;  CALL  BCL1 
IF(SWL.EQ.2)  CALL  BCL2 (PSEXIT , PREF) 
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IF(SWL.EQ.3 

IF(SWL.EQ.4 

IF(SWL.EQ.5 

IF(SWR.EQ.l 

IF(SWR.EQ.2 

IF(SWR 

IF(SWR 

IF(SWR 

IF  ((N 


EQ.3 
EQ.4 
EQ.5 


CALL  BCL3(PSINL,RINL,PREF,RREF)  RCM0111C 

CALL  BCL4(PT0TIN,RT0TIN,PREF,RREF)  RCM0112C 

CALL  BCL5  RCM0113C 

CALL  BCR1  RCM0114C 

CALL  BCR2(PSEXIT,PREF)  RCM0115C 

CALL  BCR3(PSINR,RINR,PREF,RREF)  RCM0116C 

CALL  BCR4(PTOTIN,RTOTIN,PREF,RREF)  RCM0117C 

CALL  BCR5  RCM0118C 

EQ. (50*QPRINT)).AND. (N.GE.O))  CALL  PLOT2(N,K)  RCM0119C 

40  CONTINUE  RCM0120C 

CALL  ENDPL(O)  RCM0121C 

CALL  DONEPL  RCM0122C 

STOP  RCM0123C 

END  RCM0124C 

SUBROUTINE  GLIMM(QSTOP , PSTAR, USTAR,ASTAR)  RCM0125C 

INTEGER  Q,QSTOP  RCMO 12 60 

REAL  ML , MR , MLN , MRN  '  RCM0127C 

COMMON / GLIMM1 / PGLIM , RGLIM , UGLIM , PL , RL , UL , PR , RR , UR , AL , AR , GL , GR , EPS  RCMO 12 8 0 

COMMON /GLIMM2/DT,DX, XI  RCM0129C 

DATA  Q, ML, MR/0, 100. ,100./  RCM0130C 

PSTAR=0.5*(PL+PR)  RCMO 13 10 

COEFL=SQRT(PL*RL)  RCMO 13 20 

COEFR=SQRT_(PR*RR)  RCMO  13  30 

ALPHA=1.  RCM0134G 

BEGIN  GODUNOV  ITERATION  RCM0135C 

30  Q=Q+1  RCM0136C 

IF(PSTAR.LT.EPS)  PSTAR=EPS  RCM0137C 

COMPUTE  NEXT  ITERATION  FOR  ML  AND  MR  RCMO 1380 

MLN  =  COEFL*PHI ( PSTAR , PL )  RCMO 13  9  0 

MRN=COEFR*PHI( PSTAR, PR)  RCMO 1400 

DIFML=ABS(MLN-ML)  RCM01410 

DirMR=ABS(MRN-MR)  RCM01420 

ML=MLN  RCM0143C 

MR=MRN  RCMO 1440 

COMPUTE  NEW  PSTAR  RCMO 1450 

PTIL=PSTAR  RCM0146C 

PSTAR= (UL-UR+PL/ML+PR/MR)/ (1. /ML+1. /MR)  RCMO 1470 

PSTAR=ALPHA*PSTAR+(1.-ALPHA)*PTIL  RCM01480 

IF(Q.LE.QSTOP)  GOTO  10  RCM01490 

IF(ABS(PSTAR-PTIL) .LT.EPS)  GOTO  20  RCM01500 

COMPUTE  NEW  ALPHA  RCMO 15 10 

ALPHA =0.5 -ALPHA  RCM01520 

Q=0  RCM01530 

IF((1. -ALPHA) .LT.EPS)  GOTO  20  RCM0154C 

10  IF (DIFML.GE.EPS)  GOTO  30  RCMO1550 

IF(DIFMR.GE.EPS)  GOTO  30  RCMO1560 

END  OF  GODUNOV  ITERATION;  COMPUTE  USTAR  RCMO 1570 

20  USTAR=(PL-PR+ML*UL+MR*UR)/(ML+MR)  RCMO 1580 

BEGIN  SAMPLING  PROCEDURE  RCM0159C 

IF  (XI.LT.USTAR-DT)  GO  TO  40  RCMO 1600 

RIGHT  SIDE;  SELECT  CASE  OF  SHOCK  OR  EXPANSION  RCMO 16 10 

IF  (PSTAR. LT. PR)  GO  TO  50  RCM0162C 

RIGHT  WAVE  IS  A  SHOCK  WAVE  RCM0163C 

WR=UR+MR/RR  RCM0164C 
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W: 


60 


50 


70 


80 


40 


100 


90 


IF  (XI.LT.WR*DT)  GO  TO  60 

RIGHT  OF  RIGHT  SHOCK  CASE 

RGLIM=RR 

PGLIM=PR 

UGLIM=UR 

RETURN 

LEFT  OF  RIGHT  SHOCK  CASE 

RGLIM=-MR/ (USTAR-WR) 

PGLIM=PSTAR 

UGLIM=USTAR 

RETURN 

RIGHT  WAVE  IS  A  RAREFACTION  WAVE 

CONST=PR/RR**GR 

RSTAR= ( PSTAR/ CONST ) ** ( 1 . / GR ) 

ASTAR=  SQRT ( GR*PSTAR/ RSTAR ) 

AR=SQRT(GR*PR/RR) 

IF  (XI.GE. (USTAR+ASTAR)*DT)  GO  TO  70 

LEFT  OF  RIGHT  FAN  CASE 

RGLIM= RSTAR 

UGLIM=USTAR 

PGLIM= PSTAR 

RETURN 

SELECT  RIGHT  OF  FAN 

IF  (XI.GE. (UR+AR)*DT) 

IN  RIGHT  FAN  CASE 

UGLIM=2. / (GR+1. )*(XI/DT 

RGLIM=(((AR+(GR-1. )/2.* 

PGLIM=  CONST  *RGLIM**GR 

RETURN 

RIGHT  OF  RIGHT  FAN  CASE 

RGLIM=RR 

PGLIM=PR 

UGLIM=UR 

RETURN 

LEFT  SIDE 

IF  (PSTAR 

LEFT  WAVE 

WL=UL-ML/RL 

IF  (XI.GE.WL-DT)  GO  TO  100 

LEFT  OF  LEFT  SHOCK  CASE 

RGLIM=RL 

PGLIM=PL 

UGLIM=UL 

RETURN 

RIGHT  OF  LEFT  SHOCK  CASE 

RGLIM=ML/ (USTAR-WL) 

PGLIM=PSTAR 

UGLIM=USTAR 

RETURN 

LEFT  WAVE  IS  A  RAREFACTION  WAVE 

CONST=PL/RL**GL 

RSTAR=( PSTAR/ CONST )**(1./GL) 

ASTAR=  SQRT ( GL- PSTAR/ RSTAR ) 

AL=SQRT(GL-VPL/RL) 


OR  IN  FAN 
GO  TO  80 


-AR+(GR-1. )/2." 
(UGLIM-UR))**2 


:UR) 
)/(GR*CONST))" 


(l./(GR-l.)) 


SELECT  CASE  OF 
LT.PL)  GO  TO  90 
IS  A  SHOCK  WAVE 


SHOCK  OR  RAREFACTION 
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GO  TO  110 


FAN  CASE 
TO  120 

1. )/2.*UL+XI/DT) 


IF  (XI.LT. (USTAR-ASTAR)*DT) 

RIGHT  OF  LEFT  FAN  CASE 

RGLIM=RSTAR 

PGLIM=PSTAR 

UGLIM=USTAR 

RETURN 

SELECT  LEFT  OF  FAN  OR  IN 
110  IF  (XI.LT.  (UL-AL)-'DT)  GO 

IN  LEFT  FAN  CASE 

UGLIM=2./(GL+1. )*(AL+(GL- 

RGLIM= ( ( (AL+ (GL- 1 . ) / 2 . * (UL-UGLIM) )**2 . ) / (GL-CONST ) )** ( 1 . / (GL- 1 . ) ) 

PGLIM=CONST*RGLIM**GL 

RETURN 

LEFT  OF  LEFT  FAN  CASE 
120  RGLIM=RL 

PGLIM=PL 

UGLIM=UL 

RETURN 

END 

FUNCTION  PHI(Y,Z) 

REAL  RMU 

COMMON / FUN 1 / G , PA , RA , UA , RB , RMU 

EPS=l.E-06- 

PARAM=Y/Z 

IF  (ABS(l.-PARAM) .GE.EPS)  GO  TO  10 

PHI=SQRT(G) 

RETURN 
10  IF  (PARAM.GE.l. )  GO  TO  20 

PHI= (G-l . ) /2 . *(1. - PARAM) / ( SQRT(G)* ( 1 . - PARAM** ( (G-l . ) / (2 . *G ) ) ) ) 

RETURN 
20  PHI=SQRT((G+1. )/2 .*PARAM+ (G-l. )/2. ) 

RETURN 

END 

FUNCTION  PHIl(PB) 

REAL  RMU 

COMMON/ FUN1/G , PA , RA , UA , RB , RMU 

PHI1=  (PB-PA)*SQRT( (1. -RMU**2.  )  /  (RA*  (PB  +  RMU**2  .  *PA)  ) ) 

RETURN 

END 

FUNCTION  PSI(PB) 

REAL  RMU 

COMMON/ FUN1 / G , PA , RA , UA , RB , RMU 

PSI= SQRT ( 1 . -RMU— 4 . ) /RMU**2 . / SQRT (RA)*PA** ( 1 . / (2 . *G) )*(PB** ( (G- 1 . 
*/(2.*G))-PA**((G-l.)/(2.*G))) 

RETURN 

END 

SUBROUTINE  INITl 

DIMENSION  P(203),R(203),U(203),A(203),S(203),X(203) 

COMMON/ FUN1/ G , PA , RA , UA , RB , RMU 

COMMON/ SUB  S/P,R,U,A,S,X 

DO  10  1=1,9,2 

P(I)=810600.00 

R(I)=0.7132 

U(I)=644.4 


RCM0219C 
RCM0220C 
RCM0221C 
RCM0222C 
RCM02230 
RCM02240 
RCM02250 
RCM02260 
RCM02270 
RCM02280 
RCM02290 
RCM02300 
RCM02310 
RCM02320 
RCM02330 
RCM02340 
RCM02350 
RCM02360 
RCM02370 
RCM02380 
RCM02390 
RCM02400 
RCM02410 
RCM02420 
RCM0243G 
RCM02440 
RCM02450 
RCM0246G 
RCM02470 
RCM02480 
RCM02490 
RCM02500 
RCM02510 
RCM02520 
RCM02530 
RCM02540 
RCM02550 
RCM02560 
RCM02570 
RCM02580 
RCM02590 
RCM02600 
)RCM02610 
RCM02620 
RCM02630 
RCM02640 
RCM02650 
RCM02660 
RCM02670 
RCM02680 
RCM02690 
RCM02700 
RCM02710 
RCM02720 
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A(I)=SQRT(G*P(I)/R(I)) 

10  CONTINUE 

DO  20  1=11,203,2 

P(I)=101325.0 

R(I)=1.22 

U(I)=0.0 

A(I)=SQRT(G*P(I)/R(I)) 
20  CONTINUE 

RETURN 

END 

SUBROUTINE  INIT2R(PSEXIT , PREF , RREF ) 

DIMENSION  P(203),R(203),U(203),A(203),S(203),X(203) 

COMMON/ FUN1/ G , PA , RA , UA , RB , RMU 

COMMON / SUB  S/P,R,U,A,S,X 

DO  10  1=3,201,2 

P(I)=PREF 

R(I)=RREF 

U(I)=0.0 

A(I)=SQRT(P(I)*G/R(I)) 
10  CONTINUE 

P(D  =  P(3) 

R(1)=R(3) 

U(l)=-U(3)  . 

A(1)=SQRT(G*P(1)/R(1)) 

P(203)=PSEXIT 

R(203)=R(201) 

PA=P(201) 

RA=R(201) 

UA=U(201) 

PB=P(203) 

RB=R(203) 

IF(PA.GT.PB)  GO  TO  20 

U(203)=UA-PHI1(PB) 

GO  TO  30 
20  U(203)=UA-PSI(PB) 
30  A(203)=SQRT(G*P(203)/R(203)) 

RETURN 

END 

SUBROUTINE  INIT2L(PSEXIT) 

DIMENSION  P(203),R(203),U(203),A(203),S(203),X(203) 

COMMON / FUN1 / G , PA , RA , UA , RB , RMU 

COMMON / SUBS / P , R , U , A , S , X 

DO  10  1=3,201,2 

P(I)=285080.0 

R(I)=0.897 

U(I)=0.0 

A(I)=SQRT(G-P(I)/R(I)) 
10  CONTINUE 

RETURN 

END 

SUBROUTINE  INIT3L(PSINL,RINL) 

DIMENSION  P(203),R(203),U(203),A(203),S(203),X(203) 

COMMON/ FUN1/ G , PA , RA , UA , RB , RMU 

COMMON / SUBS / P , R , U , A , S , X 
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DO  10  1=3,201,2 

P(I)=2390000.0 

R(I)=9.787 

U(I)=0.0 

A(I)=SQRT(G*P(I)/R(I)) 
10  CONTINUE 

P(1)=PSINL 

R(1)=RINL 

PA=P(3) 

RA=R(3) 

UA=U(3) 

PB=P(1) 

U(1)=UA+PHI1(PB) 

A(1)=SQRT(G*P(1)/R(1)) 

P(203)=P(201) 

R(203)=R(201) 

U(203)=-U(201) 

A( 203 )= SORT (G*P ( 203 )/R( 203)) 

RETURN 

END 

SUBROUTINE  INIT3R(PSINR,RINR) 

DIMENSION  P(203),R(203),U(203),A(203),S(203),X(203) 

COMMON/ SUBS/ P,R,U, A, S,X 

COMMON/ FUN1/ G , PA , RA , UA , RB , RMU 

DO  10  1=3,201,2 

P(I)=2421667.5 

R(I)=9.787 

U(I)=0.0 

A(I)=SQRT(G*P(I)/R(I)) 
10  CONTINUE 

P(203)=PSINR 

PA=P(201) 

RA=R(201) 

UA=U(201) 

PB=P(203) 

U(203)=UA-PHI1(PB) 

R(203)=RINR 

A(203)=SQRT(G*P(203)/R(203)) 

P(1)=P(3) 

R(1)=R(3) 

U(l)=-U(3) 

A(1)=SQRT(G*P(1)/R(1)) 

RETURN 

END 

SUBROUTINE  BCL1 

DIMENSION  P(203),R(203),U(203),A(203),S(203),X(203) 

COMMON / SUB  S/P,R,U,A,S,X 

COMMON/ FUN1/ G , PA , RA , UA , RB , RMU 

P(1)=P(3) 

R(1)=R(3) 

U(l)=-U(3) 

A(1)=SQRT(G-P(1)/R(1)) 

RETURN 

END 


RCM03270 
RCM03280 
RCM03290 
RCM03300 
RCM03310 
RCM03320 
RCM03330 
RCM03340 
RCM03350 
RCM03360 
RCM03370 
RCM03380 
RCM03390 
RCM03400 
RCM03410 
RCM03420 
RCM03430 
RCM03440 
RCM03450 
RCM03460 
RCM03470 
RCM03480 
RCM03490 
RCM03500 
RCM03510 
RCM03520 
RCM0  3  5  30 
RCM03540 
RCM03550 
RCM03560 
RCM03570 
RCM03580 
RCM03590 
RCM03600 
RCM03610 
RCM03620 
RCM03630 
RCM03640 
RCM03650 
RCM03660 
RCM03670 
RCM03680 
RCM03690 
RCM03700 
RCM03710 
RCM03720 
RCM03730 
RCM03740 
RCM03750 
RCM03760 
RCM03770 
RCM03780 
RCM03790 
RCM03800 
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\V 


10 
20 


SUBROUTINE  BCRl 

DIMENSION  P(203),R(203),U(203) , A(203 ) , S (203 ) ,X(203 ) 

COMMON/ SUBS /P,R,U, A, S,X 

COMMON/ FUN1/ G , PA , RA , UA , RB , RMU 

P(203)=P(201) 

R(203)=R(201) 

U(203)=-U(201) 

A(203)=SQRT(G*P(203)/R(203)) 

RETURN 

END 

SUBROUTINE  BCL2 (PSEXIT , PREF) 

DIMENSION  P(203),R(203),U(203),A(203),S(203),X(203) 

COMMON/ SUBS /P,R,U, A, S,X 

COMMON/ FUN1/ G , PA , RA , UA , RB , RMU 

P(1)=PSEXIT/PREF 

R(1)=R(3) 

U(1)=U(3)« 

A(1)=SQRT(G«P(1)/R(1)) 

RETURN 

END 

SUBROUTINE  BCR2 (PSEXIT , PREF) 

DIMENSION  P(203),R(203),U(203),A(203),S(203),X(203) 

COMMON/ SUBS /P,R,U, A, S,X 

COMMON/ FUNl/G , PA, RA , UA , RB , RMU 

P(203)=PSEXIT/PREF 

R(203)=R(201) 

PA=P(20l) 

RA=R(201) 

UA=U(201) 

PB=P(203) 

RB=R(203) 

IF(PA.GT.PB)  GO  TO  10 

U(203)=UA-PHI1(PB) 

GO  TO  20 

U(203)=UA-PSI(PB) 

A(203)=SQRT(G*P(203)/R(203)) 

RETURN 

END 

SUBROUTINE  BCL3 ( PSINL , RINL , PREF , RREF ) 

DIMENSION  P(203),R(203),U(203),A(203),S(203),X(203) 

COMMON / SUBS / P , R , U , A , S , X 

COMMON/ FUN1/ G , PA , RA , UA , RB , RMU 

P(1)=PSINL/PREF 

R(1)=RINL/RREF 

PA=P(3) 

RA=R(3) 

UA=U(3) 

PB=P(1) 

U(1)=UA+PHI1(PB) 

A( 1 ) = SQRT (G*P ( 1 ) /R( 1 ) ) 

RETURN 

END 

SUBROUTINE  BCR3 ( PSINR , RINR , PREF , RREF ) 

DIMENSION  P(203),R(203),U(203)  ,A'(203)  ,  S(203)  ,X(203) 
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60 


50 


COMMON/ SUBS /P,R,U, A, S,X 

COMMON/ FUN1 / G , PA , RA , UA , RB , RMU 

P(203)=PSINR/PREF 

R(203)=RINR/RREF 

PA=P(201) 

RA=R(201) 

UA=U(201) 

PB=P(203) 

U(203)=UA-PHI1(PB) 

A(203)=SQRT(G*P(203)/R(203)) 

RETURN 

END 

SUBROUTINE  BCL4  (PTOTIN,RTOTIN , PREF ,RREF) 

INTEGER  QOUT 

DIMENSION  P(203),R(203),U(203),A(203),S(203),X(203) 

DIMENSION  XARRAY(IOO) 

COMMON / SUB  S/P,R,U,A,S,X 

COMMON / FUN 1 / G , PA , RA , UA , RB , RMU 

COMMON  XARRAY,N1 

N1=N1+1 

QOUT=Nl/5 

PTOT=PTOTIN/PREF 

RTOT=RTOTIN/RREF 

ATOT=SQRT(G*PTOT/RTOT) 

STOT=ALOG(PTOT/RTOT**G) 

U(1)=U(3) 

A ( 1 )  =  SQRT ( ATOT * * 2 . - ( G -  1 . ) / 2 . * ABS (U ( 1 ) ) **2 . ) 

AMACH=U(1)/A(1) 

IF(AMACH.LT.O.O)  GO  TO  60 

P(l)=PTOT/(l.+(G-l. ) / 2 . * AMACH** 2 . )**(G/(G-1. )) 

R(l)=RTOT/(l.+(G-l. )  / 2 . *AMACH**2 . )**(1./(G-1. )) 

S ( 1 ) = ALOG ( P ( 1 ) / R ( 1 ) **G ) 

GO  TO  50 

P(l)=PTOT/(l.+(G-l. )/2.*ABS(AMACH)**2. ) 

R(l)=RTOT/ (l.+(G-l. )/2.*ABS(AMACH)**2. ) 

S(l)=ALOG(P(l)/R(l)**G) 

RETURN 

END 

SUBROUTINE  BCR4  (PTOTIN ,RTOTIN , PREF ,RREF) 

INTEGER  QOUT 

DIMENSION  P(203),R(203),U(203),A(203),S(203),X(203) 

DIMENSION  XARRAY(IOO) 

COMMON/ SUBS /P,R,U, A, S,X 

COMMON/ FUN1/ G , PA , RA ,UA , RB , RMU 

COMMON  XARRAY,N1 

N1=N1+1 

QOUT=Nl/25 

PTOT=PTOTIN/PREF 

RTOT=RTOTIN/RREF 

ATOT=SQRT(G-PTOT/RTOT) 

STOT  =  ALOG(PTOT/RTOT'v"G) 

U(203)=U(201) 

A(  203 )  =  SQRT (ATOT **2 . - (G 

AMACH=U(203)/A(203) 


(G/(G-1.)) 
(l./(G-l.)) 


l.)/2.*ABS(U(203))**2.) 


RCM04350 
RCM04360 
RCM04370 
RCM04380 
RCM04390 
RCM04400 
RCM04410 
RCM04420 
RCM04430 
RCM04440 
RCM04450 
RCM04460 
RCM04470 
RCM04480 
RCM04490 
RCM04500 
RCM04510 
RCM04520 
RCM04530 
RCM04540 
RCM04550 
RCM04560 
RCM04570 
RCM04580 
RCM04590 
RCM04600 
RCM04610 
RCM04620 
RCM04630 
RCM04640 
RCM04650 
RCM04660 
RCM04670 
RCM04680 
RCM04690 
RCM04700 
RCM04710 
RCM04720 
RCM04730 
RCM04740 
RCM04750 
RCM04760 
RCM04770 
RCM04780 
RCM04790 
RCM04800 
RCM04810 
RCM04820 
RCM048  30 
RCM04840 
RCM04850 
RCM04860 
RCM048  70 
RCM04880 
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•u 


60 


50 


20 


10 


40 
50 


60 


)**(G/(G-1.)) 
)**(1./(G-1.)) 


IF(AMACH.LT.O.O)  GO  TO  60 

P(203)=PTOT/(l.+(G-l. )/2.*AMACH**2. )**(G/(G-1. )) 

R(203)=RTOT/(1.+(G-1. ) / 2 . *AMACH**2 . )**(!. /(G-l. )) 

S(203)=ALOG(P(203)/R(203)**G) 

GO  TO  50 

P(203)=PTOT/(1.+(G-1. )/2.*ABS(AMACH)**2 

R(203 ) =RTOT/ (1. + (G-l. ) /2 . *ABS (AMACH)**2 

S(203)=ALOG(P(203)/R(203)**G) 

RETURN 

END 

SUBROUTINE  BCL5 

DIMENSION  P(203),R(203),U(203),A(203),S(203),X(203) 

COMMON/ SUBS/P,R,U, A, S,X 

P(1)=P(3) 

R(1)=R(3) 

U(1)=U(3) 

A(1)=A(3) 

RETURN 

END 

SUBROUTINE  BCR5 

DIMENSION  P(203),R(203),U(203),A(203),S(203) ,X(203) 

COMMON / SUB  S/P,R,U,A,S,X 

P(203)=P(201) 

R(203)=R(201) 

U(203)=U(201) 

A(203)=A(201) 

RETURN 

END 

FUNCTION  WDP(II) 

DIMENSION  WNORM(12),IDIGT(12) 

COMMON/ S AMPLE /WNORM,IDIGT 

IF  (II.EQ.O)  GO  TO  10 

Ll  =  2 

L2=l 

DO  20  JJ=1,12 

IDIGT(JJ)=0 

WNORM( J J ) = 1 . /FLOAT (LI ** J J ) 

CONTINUE 

WDP=0. 

RETURN 

DO  40  JJ=1,12 

Ll=2 

L2=l 

KJO=IDIGT(JJ) 

KJN=MOD((KJ0+l) ,L1) 

IDIGT(JJ)=KJN 

IF  (KJO.LT.KJN)  GO  TO  50 

CONTINUE 

SUM=0. 

DO  60  JJ=1,12 

KNEW=MOD ( IDIGT ( J J ) *L2 , LI ) 

SUM=  SUM+  FLOAT ( KNEW ) * WNORM ( J J ) 

CONTINUE 

WDP=SUM 
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•it.' 


RETURN 

END 

SUBROUTINE  PLOTl(K) 

DIMENSION  XORG(4) ,Y0RG(4) , YMAX(4) , YMIN (4) 


5,4.75,0.5,4.75/ 
5,0.5,4.75,4.75/ 
50,3.0,0.5,2.0/ 
5,0.0,-0.5,-2.0/ 


DATA  XORG/0 
DATA  YORG/0 
DATA  YMAX/3 
DATA  YMIN/0 
DO  10  1=1,4 

CALL  PHYSOR(XORG(I) ,YORG(I)) 
AREA2D(3.5,3.5) 


CALL 
CALL 
CALL 
CALL 


FRAME 
GRAF(0. , 
ENDGR(O) 
10  CONTINUE 

CALL  PHYSOR(8 
AREA2D(2 
FRAME 
GRAF(0. , 
ENDGR(O) 


SCALE ',1.0, YMIN ( I ) , ' SCALE ' , YMAX ( I ) ) 


5,0.5) 
25,7.75) 


SCALE' ,1. ,0, 'SCALE' ,K) 


CALL 
CALL 
CALL 
CALL 
RETURN 
END 

SUBROUTINE- PLOT2(N,K) 

DIMENSION  XORG(4) ,Y0RG(4) ,YMAX(4) ,YMIN(4) ,KNT(4) ,IYNAM(10) 
DIMENSION  PARRAY ( 100 ) , RARRAY ( 100 ) , UARRAY ( 100 ) , SARRAY ( 100 ) , XARRAY ( 
»00 ) 
DIMENSION  P(203) ,R(203),U(203),A(203) ,S(203) ,X(203) 
COMMON/ SUBS /P,R,U, A, S,X 
COMMON  XARRAY 
DATA  XORG/0. 5, 4. 75, 0.5, 4 

YORG/0. 5, 0.5, 4. 75, 4 

YMAX/ 3. 50, 3. 0,0.  5,  2 

YMIN/ 0.5, 0.0,  -0.5 

KNT/1,4,6,9/ 

IYNAM/ 'PRES' 
*, 'ENTR' , 'OPY$'/ 
DO  200  1=1,100 
PARRAY(I)=P (1-2+1) 
RARRAY(I)=R(I*2+1) 
UARRAY(I)=U(I*2+1) 
SARRAY(I)=S (1*2+1) 
200  CONTINUE 

DO  300  1=1,4 

CALL  PHYSOR(XORG(I) ,YORG(I)) 

AREA2D(3.5,3.5) 

XNAME ( ' X ' , 1 ) 

YNAME ( IYNAM (KNT ( I ) ) , 100 ) 

GRAF(0. , 'SCALE' ,1.0,YMIN(I) 

EQ.l)  CALL  SETCLR( 'YELLOW ) 


DATA 
DATA 
DATA 
DATA 
DATA 


75/ 
75/ 
0/ 
2.0/ 


SURE ' , ' $ 


DENS ' , ' ITY$ ' , ' VELO ' , ' CITY ' , ' $ 


CALL 
CALL 
CALL 
CALL 
IF(I 
IF(I 
IF(I 
IF(I 
IF(N 
IF(I 


SCALE' ,YMAX(I)) 


EQ.2)  CALL  SETCLR( 'CYAN' ) 

EQ.3)  CALL  SETCLR( 'RED' ) 

EQ.4)  CALL  SETCLR( 'MAGENTA' ) 

EQ . K )  CALL  SETCLR ( ' WHITE ' ) 

EQ.l)  CALL  CURVE  ( XARRAY, PARRAY, 100 , 0 ) 


RCM0543C 
RCM05440 
RCM05450 
RCM05460 
RCM05470 
RCM05480 
RCM05490 
RCM05500 
RCM05510 
RCM05520 
RCM05530 
RCM05540 
RCM05550 
RCM05560 
RCM05570 
RCM05580 
RCM05590 
RCM05600 
RCM05610 
RCM05620 
RCM05630 
RCM05640 
RCM05  6  50 
RCM05660 

1RCM05670 
RCM05680 
RCM05690 
RCM05700 
RCM05710 
RCM05720 
RCM05730 
RCM05740 
RCM05750 
RCM05760 

'RCM05770 
RCM05780 
RCM05790 
RCM05800 
RCM05810 
RCM05820 
RCM05830 
RCM05840 
RCM05850 
RCM05860 
RCM05870 
RCM05880 
RCM05890 
RCM05900 
RCM05910 
RCM05920 
RCM05930 
RCM05940 
RCM05950 
RCM05960 
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IF(I.EQ.2)  CALL  CURVE  ( XARRAY , RARRAY , 100 , 0 ) 

IF(I.EQ.3)  CALL  CURVE  ( XARRAY, UARRAY , 100 , 0 ) 

IF(I.EQ.4)  CALL  CURVE  (XARRAY , SARRAY , 100 , 0 ) 

CALL  ENDGR(O) 
300  CONTINUE 

RETURN 

END 

SUBROUTINE  GE ( SWL , SWR , N , TTOTAL , TIME , UEXMAX , PTOTIN , PREF ) 

INTEGER  SWL, SWR 

DIMENSION  P(203),R(203) ,U(203 ) ,A(203 ) , S (203 ) ,X(203 ) 

COMMON / SUBS / P , R , U , A , S , X 
C* -CALCULATION  STARTS  AT  EXHAUST  PORT  OPENING.  SUBROUTINE  STRUCTURED 
C* -ACCORDINGLY. 

IF((SWL.EQ.1).AND. (SWR.EQ.2))  GO  TO  10 

IF((SWL.EQ.4) .AND. (SWR.EQ.2))  GO  TO  30 

IF((SWL.EQ.4).AND. (SWR.EQ.l))  GO  TO  50 

IF ( ( SWL . EQ . 1 ) . AND . ( SWR . EQ . 1 ) )  RETURN 
10  PWALL=P(2) 

(PTOTIN/PREF))  GO  TO  20 


LE 


IF(PWALL 

RETURN 
20  SWL=4 

WRITE (6, 74) 

WRITE (6, 75 J  N, TTOTAL, TIME 

RETURN 
30  UEXIT=U(202) 

IF (UEXMAX. LT.UEXIT)  UEXMAX=UEXIT 

IF(UEXIT.LT.UEXMAX/2. )  GO  TO  40 

RETURN 
40  SWR=1 

WRITE (6, 76) 

WRITE (6, 75)  N, TTOTAL, TIME 

RETURN 
50  P1SH0K=P(2) 

IF(P1SH0K.GT. PTOTIN/PREF)  GO  TO  60 

RETURN 
60  SWL=1 

WRITE (6, 77) 

WRITE (6, 75)  N, TTOTAL, TIME 

74  FORMAT ( 5 X, 'INLET  PORT  OPENS  AT:') 

75  F0RMAT(5X,I4,5X,2F14.7) 

76  FORMAT ( 5 X, 'EXHAUST  PORT  CLOSES  AT:') 
7  7  FORMAT ( 5 X, 'INLET  PORT  CLOSES  AT:') 

RETURN 
END 

SUBROUTINE  SPCTRA  (N ,  SWL ,  SWR ,  TIME  ,  UEXMAX ,  PSEXIT  ,  PSOUT1 ,  PSOUT2  ,  PSOUTRCMCji 
*3  , JCOUNT , QPRINT , TTOTAL , KCOUNT ) 
INTEGER  SWL, SWR, QPRINT 

DIMENSION  P(203) ,R(203),U(203),A(203),S(203),X(203) 
COMMON/ SUBS /P,R,U, A, S,X 
C-CALCULATION  STARTS  AT  HP  GAS  IN  PORT.  JCOUNT 
IF((SWL.EQ.1).AND. (SWR.EQ.3))  GO  TO  10 
IF((SWL.EQ.2) .AND. (SWR.EQ.3))  GO  TO  20 
IF( (SWL. EQ. 2) .AND. (SWR.EQ.l))  GO  TO  30 
IF( ( SWL. EQ. 5) .AND. (SWR.EQ.l))  GO  TO  40 


IS  NUMBERED  ACCORDINGLY- 
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r;..l 


50 
60 
70 


10 


11 


12 
13 

20 


22 


23 


24 
25 

26 

30 


31 


32 
33 

40 


41 


IF((SWL.EQ.2).AND. (SWR.EQ.5))  GO  TO 

IF( (SWL. EQ. 2). AND. (SWR.EQ.4))  GO  TO 

IF((SWL.EQ. 1) .AND. (SWR.EQ.4))  GO  TO 

IF(U(3) .LT.0.0)  GO  TO  11 

RETURN 

JC0UNT=JC0UNT+1 

PSEXIT=PS0UT1 

SWL=2 

WRITE(6,12) 

WRITE ( 6 , 13 )N , TIME , SWL , SWR , JCOUNT 

FORMAT (5X, 'H.P.  AIR  OUT  PORT  OPENS  AT') 

F0RMAT(5X,I4,5X,F9.7,5X,3I3) 

RETURN 

DO  26  1=5,199,2 

IF((R(I)-R(I+2)).GT.0.1)  GO  TO  22 

GO  TO  26 

XCNTCT=X(I) 

UCNTCT=U(I) 

TCNTCT=XCNTCT/ ABS (UCNTCT ) 

AHEAD=A(199) 

THEAD=1.0/A(199) 

IF(TCNTCT.LE.THEAD)  GO  TO  23 

RETURN 

JCOUNT=JCOUNT+l 

SWR=1 

WRITE (6, 24) 

WRITE ( 6 , 25 ) N , TIME , SWL , SWR , JCOUNT 

FORMAT (5X, 'H.P.  GAS  IN  PORT  CLOSES  AT') 

FORMAT ( 5X , 14 , 5X , F9 . 7 , 5X , 31 3 ) 

RETURN 

CONTINUE 

IF(JC0UNT.EQ.4) 

IF(JC0UNT.EQ.6) 

IF(JC0UNT.EQ.8) 

IF((R(2)-R(4)) 

RETURN 

JCOUNT=JCOUNT+l 

SWL=5 

WRITE(6,32) 

WRITE ( 6  ,  3  3  )  N , TIME , SWL , SWR , JCOUNT 

FORMAT (5X, 'HP  AIR  OUT  PORT  CLOSES 

FORMAT ( 5X , 14 , 5X , F9 . 7 , 5X , 3 13 ) 

RETURN 

GO 

GO 


GO  TO  80 
GO  TO  90 
GO  TO  100 
GT.0.1)  GO  TO  31 


AND  TUNING  PORT  LI  OPENS  AT') 


7) 
0) 


TO 
TO 


110 
41 


42 
43 


IF(JCOUNT.EQ 

IF(U(3).GE.O 

RETURN 

JCOUNT=JCOUNT+l 

PSEXIT=PS0UT2 

SWL=2 

WRITE (6, 42) 

WRITE (6, 43 )N, TIME, SWL, SWR, JCOUNT 

FORMAT (5X, 'TUNING  PORT  LI  CLOSES 

FORMAT ( 5X , 14 , 5X , F9 . 7 , 5X , 313 ) 

RETURN 


AND  EXHAUST  PORT  El  OPENS  AT') 


RCM0651C 
RCM0652C 
RCM0653C 
RCM0654C 
RCM0655C 
RCM0656C 
RCM0657C 
RCM0658C 
RCM06  5  9C 
RCM0660C 
RCM06  61C 
RCM06620 
RCM0663C 
RCM06640 
RCM06650 
RCM06660 
RCM06670 
RCM0668C 
RCM06690 
RCM06700 
RCM06710 
RCM06720 
RCM06730 
RCM0674C 
RCM06750 
RCM06760 
RCM06770 
RCM06780 
RCM06790 
RCM06800 
RCM06810 
RCM06820 
RCM068  30 
RCM06840 
RCM06850 
RCM06860 
RCM06870 
RCM06880 
RCM06890 
RCM06900 
RCM06910 
RCM06920 
RCM06930 
RCM06940 
RCM06950 
RCM06960 
RCM06970 
RCM06980 
RCM06990 
RCM07000 
RCM07010 
RCM07020 
RCM07030 
RCM07040 
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,  if 


80  IF(ABS(U(201)).GT. .0001)  GO  TO  81 
RETURN 

81  JC0UNT=JC0UNT+1 
SWR=5 

WRITE (6, 82) 

WRITE ( 6 , 8  3 ) N , TIME , SWL , SWR , JCOUNT 

82  FORMAT (5X, 'TUNING  PORT  Rl  OPENS  AT') 

83  F0RMAT(5X,I4,5X,F9.7,5X,3I3) 
RETURN 

50  IF(JCOUNT.EQ.9)  GO  TO  120 
IF(ABS(U(201)-U(3)).LE. 0.001)  GO  TO  51 
RETURN 

51  JCOUNT=JCOUNT+l 
SWR=1 

WRITE (6, 52) 

WRITE ( 6 , 5  3 ) N , TIME , SWL , SWR , JCOUNT 

52  FORMAT ( 5 X, 'TUNING  PORT  Rl  CLOSES  AT') 

53  F0RMAT(5X,I4,5X,F9.7,5X,3I3) 
RETURN 

90  THEAD1=X(201)/(ABS(U(3))+A(3)) 
KC0UNT=KC0UNT+1 
IF(KCOUNT.EQ.l)  TT0T1=TT0TAL 
IF(TTOTAL.GE.  (TTOTl+THEADl) )  GO  TO  91 
RETURN 

91  JCOUNT=JCOUNT+l 
SWL=5 
WRITE(6,92) 

WRITE ( 6 , 9  3 ) N , TIME , SWL , SWR , JCOUNT 

92  FORMAT (5X, 'EXHAUST  PORT  El  CLOSES  AND  TUNING  PORT  L2  OPENS  AT') 

93  F0RMAT(5X,I4,5X,F9.7,5X,3I3) 
RETURN 

110  IF(U(3).GE.0.0)  GO  TO  111 
RETURN 

111  JC0UNT=JC0UNT+1 
PSEXIT=PS0UT3 
SWL=2 

WRITE (6, 112) 

WRITE ( 6 , 1 13 ) N , TIME , SWL , SWR , JCOUNT 

112  FORMAT (5X, 'TUNING  PORT  L2  CLOSES  AND  EXHAUST  PORT  E2  OPENS  AT') 

113  F0RMAT(5X,I4,5X,F9.7,5X,3I3) 
RETURN 

100  IF(ABS(U(201)).GT. .0001)  GO  TO  101 
RETURN 

101  JCOUNT=JCOUNT+l 
SWR=5 
WRITE(6,102) 

WRITE ( 6 , 10  3 )N , TIME , SWL , SWR , JCOUNT 

102  FORMAT (5X, 'TUNING  PORT  R2  OPENS  AT') 

103  F0RMAT(5X,I4,5X,F9.7,5X,3I3) 
RETURN 

120  IF(ABS(U(201)-U(3)).LE. 0.0001)  GO  TO  121 
RETURN 

121  JC0UNT=JC0UNT+1 
SWR=4 
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WRITE(6,122) 

WRITE ( 6  ,  12  3  )  N , TIME , SWL , SWR , JCOUNT 

122  FORMAT ( 5 X, 'TUNING  PORT  R2  CLOSES  AND  L.P 

123  F0RMAT(5X,I4,5X,F9.7,5X,3I3) 
RETURN 

60  IF((R(4)-R(2)).GT.0.1)  GO  TO  61 
RETURN 

61  JCOUNT=JCOUNT+l 
SWL=1 

WRITE (6, 62) 

WRITE ( 6 , 6  3 ) N , TIME , SWL , SWR , JCOUNT 

62  FORMAT (5X, 'EXHAUST  PORT  E2  CLOSES  AT') 

63  F0RMAT(5X,I4,5X,F9.7,5X,3I3) 
RETURN 

70  IF(U(201) .GE.0.0)  GO  TO  71 
RETURN 

71  JC0UNT=0 
SWR=1 
WRITE (6, 72) 

WRITE (  6  ,  7  3  )  N  , TIME , SWL , SWR , JCOUNT 

72  FORMAT (5X, 'CYCLE  COMPLETED.') 

73  F0RMAT(5X,I4,5X,F9.7,5X,3I3) 
RETURN 

END 


AIR  INLET  OPENS  AT' ) 


RCM07590 
RCM07600 
RCM07610 
RCM07620 
RCM07630 
RCM07640 
RCM07650 
RCM07660 
RCM07670 
RCM07680 
RCM07690 
RCM07700 
RCM07710 
RCM07720 
RCM07730 
RCM07740 
RCM07750 
RCM07760 
RCM07770 
RCM07780 
RCM07790 
RCM07800 
RCM07810 
RCM07820 
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APPENDIX  B 
PROGRAM  RCM 
B.  1 .   Program  Description 
B.  1.1.   Computational  Grid 

The  computational  region  is  divided  into  100  cells;  the  solution  grid 
points  are  odd  numbered,  e.g.,  3,  5,  7  ...,  201  with  1  and  203  being  the 
points  where  the  boundary  conditions  are  specified.   The  even  numbered  points, 
2,  4,  6  ...,  202  are  intermediate  locations  where  solutions  are  stored  before 
being  assigned  to  the  solution  grid  points.   See  Fig.  (2). 

B.1.2.   Data  Input 

Data  for  various  ports  (exhaust,  inlet,  etc.)  is  specified  in 
dimensional  form  in  S.I.  units  (Pascal  (N/m^)  for  pressure,   kg/ra^   for 
density,   ra/s   for  velocity  etc.).   Reference  values  are  also  specified  in 
like  manner.   See  lines  RCM00210  through  00250. 

Initial  data  is  specified  through  a  call  to  an  appropriate 
subroutine,  depending  on  where  the  calculation  is  started  for  a  particular 
wave  diagram.   For  the  example  given  in  section  II  on  the  Spectra  Technology 
wave  diagram,  the  computation  is  started  at  the  point  when  the  high  pressure 
gas  inlet  port  just  opens.   The  call  for  initial  data  is  made  to  subroutine 
INIT3R,  which  prescribes  data  consistent  with  a  solid  wall  boundary  at  Che 
left  and  a  'piston'  inflow  boundary  at  the  right. 

B.1.3.   Non -d i mens ionalizat ion 

Non-dimensionalization  is  carried  out  in  lines  00540  through  00610 
with  entropy  defined  as 

S  =  Jin  (£-) 
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Note  that  velocities  are  all  referred  to  a  reference  sonic  velocity  defined 
by 

a  c    =    pref 

aref     — — 

^ref 

B.1.4.   Structure 

The  main  program  loop  starts  at  line  00630,  for  the  number  of  time 
steps  specified.   The  time  step  is  computed  according  to  the  appropriate  CFL 
condition  for  the  method,  and  a  random  number  for  the  time  step  is  generated 
by  a  call  to  the  function  subroutine  WDP . 

A  secondary  loop  to  define  the  sequence  of  local  Riemann  problems  for 
the  time  step  is  set  up  at  line  00750.   For  each  Riemann  problem  defined,  a 
call  is  made  to  subroutine  GL1MM  which  i)  solves  the  Riemann  problem,  and  ii) 
samples  the  solution  using  the  random  number  generated.   The  subroutine  then 
returns  the  sampled  solution  as  the  parameters  PGLIM,  RGLIM,  UGLIM  for  the 
pressure,  density  and  velocity  respectively.   These  solutions  are  initially 
stored  in  the  even  numbered  intermediate  locations  on  the  grid,  and  are  then 
assigned  to  either  the  left  or  the  right  solution  grid  point  depending  on 
whether  the  random  number  was  In  the  negative  or  the  positive  half  of  the 
interval  respectively. 

A  call  is  then  made  to  one  of  the  modular  subroutines  structured  for 
particular  types  of  wave  diagrams,  lines  01050-01080,  and  the  others  are 
commented  out. 

Boundary  conditions  are  invoked  after  the  call  to  the  modular 
subroutines  which  return  the  proper  values  of  the  switches  SWL  and  SWR.   The 
structure  of  the  boundary  condition  subroutines  is  described  in  section  II. 
This  sequence  completes  one  pass  through  the  main  loop  and  the  process  is 
repeated  for  the  number  of  time  steps  specified. 
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B.2.   Example  Use  of  Program  RCM 

The  program  is  set  up  in  the  following  steps: 

i)    Line  00150   -  output  device  designation.  See  B.3. 

ii)   Line  00190  -  specify  the  number  of  time  steps,   k   ,  and  the 

switches  SWL  and  SWR  consistent  with  where  the 
computation  is  to  be  started. 

iii)   Lines  00210  -  prescribe  flow  data  for  various  ports  in 

through  00250  dimensional  form.   See  list  of  variables  for 
explanation  of  variable  names. 

iv)    Lines  00490  -  invoke  the  proper  initial  data  subroutine  and 
through  00530   comment  out  the  rest.   See  list  of  subroutines 

for  explanation  of  subroutine,  function  subprogram 
names . 

v)    Line  00660  -  set  the  interval  for  number  of  time  steps  at  which 

a  plot  of  the  flow  parameters  is  required. 

vi )   Lines  01050  -  user  supplied  modular  subroutine  for  particular 

through  01080  wave  diagram  to  be  computed.   Comment  out  the  rest. 

vii)   Line  01190  -   call  to  plotting  routine  should  be  consistent  with 

interval  specified  in  line  0660. 

viii)  Lines  02650  -  identify  proper  subroutine  to  prescribe  initial 

through  03700  data  (consistent  with  iv),  and  specify  the  data  in 
the  subroutine  in  dimensional  form. 

ix)    Lines  05470  -  specify  plotting  parameters,  viz.,  origins  of 
through  05990   plots,  scales,  number  of  points  to  be  plotted, 
color  of  plots,  etc.   Facility  dependent. 

The  subroutines  PL0T1  and  PL0T2  given  in  the 
listing  are  structured  for  DISSPLA  software 
installed  in  the  facility  at  NPGS. 

x)     Lines  06040  -  user  supplied  modular  subroutine  for  wave  diagram 
through  07820   to  be  computed. 


B. 3.   Execution 

The  program  is  run  in  an  interactive  mode  and  is  invoked  through  a  call 
to  DISSPLA,  available  on  most  mainframes.   After  compiling  the  program 
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(FORTRAN  H  Extended  compiler),  the  following  command  executes  it: 

DISSPLA  filename 
If  working  at  stations  equipped  with  dual  screens,  the  command  on  line  150  can 
be  of  the  type 

CALL  TEK618   +■  Tektronix  screen 
If  working  on  a  non-graphics  terminal,  or  a  single  screen  station,  this 
should  be  changed  to 

CALL  COMPRS 
which  generates  a  'DISSPLA  METAFILE'  to  be  routed  later  to  either  a  screen  or 
a  plotter,  e.g.,  VRSTEC,  IBM79,  TEK618,  etc.   Once  generated,  the  metafile  can 
be  accessed  and  routed  by  the  command 

DISSPOP  device  designation 
These  are  facility  dependent  commands  and  should  be  modified  accordingly. 

B.4.   List  of  Important  Variables   (In  Alphabetical  Order) 

A  -  sonic  velocity 

AHEAD  -  sonic  speed  of  head  wave  of  rarefaction  fan 

AMACH  -  Mach  number 

AL  -  left  side  sonic  speed  value  for  RP 

AR  -  right  side  sonic  speed  value  for  RP 

AREF  -  reference  speed  of  sound 

ASTAR  -  speed  of  sound  in  'starred'  state  of  RP  solution 

(see  Fig.  3) 

CFLNUM  -  CFL  number  for  time  step  determination 

DT  -  time  step 

DX  -  grid  cell  width 

EPS  -  small  number  for  pressure  iteration  in  RP  solver 

G  -  ratio  of  specific  heats,   y 

IDIGT  -  see  WNORM 

II  -  argument  used  in  function  subprogram  PHI  equal  to  either 

0  or  1 

JCOUNT  -  counter 

K  -  number  of  time  steps 

KCOUNT  -  counter 

N  -  counter  for  time  steps 

Nl  -  counter 
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P  -  pressure 

PA  -  flow  parameter  describing  'a'  state  in  transition  functions 

PGLIM  -  pressure  value  returned  by  subroutine  GLIMM 

PL  -  left  side  pressure  value  for  RP 

PR  -  right  side  pressure  value  for  RP 

PREF  -  reference  pressure 

PSEXIT  -  static  pressure  at  exit  or  outlet  port 

PSINL  -  static  pressure  for  incoming  'piston'  flow  on  left  side 

PSINR  -  static  pressure  for  incoming  'piston'  flow  on  right  side 

PSOUTn  -  n  =  1,2,3  -  exit  static  pressures  for  cycles  with  more  than 

one  exhaust  port 

PSTAR  -  pressure  in  'starred'  state  of  RP  solution  (see  Fig.  3) 

PTOTIN  -  total  pressure  for  isentropic  inflow 

QPRINT  -  specification  of  interval  size  for  output 

QSTOP  -  maximum  number  of  iterations  for  solution  of  Riemann 

problem,  (RP) 

R  -  density 

RA,RB  -  flow  parameters  describing  'a'  and  'b'  states  in  transition 

functions 

RGLIM  -  density  returned  by  subroutine  GLIMM 

RINL  -  static  density  for  incoming  'piston'  flow  on  left  side 

RINR  -  static  density  for  incoming  'piston'  flow  on  right  side 

RL  -  left  side  density  for  RP 

RMU  -  function  of   Y 

RR    "  right  side  density  for  RP 

RREF  -  reference  density 

RTOTIN  -  total  density  for  isentropic  inflow 

S  -  entropy 

SWL  -  switch  for  left  boundary 

SWR  -  switch  for  right  boundary 

TCNTCT  -  time  taken  by  contact  surface  to  travel  a  certain  distance 

THEAD  -  time  taken  by  head  wave  of  expansion  to  travel  a  certain 

distance 

TIME  -  real  time  in  seconds 

TIMEREF  -  reference  time 

TTOTAL  -  cumulative  non-dimensional  time  for  number  of  time  steps 

U  -  velocity 

UA  -  flow  parameter  for  'a'  state  in  transition  functions 

UCNTCT  -  velocity  of  contact  surface 

UEXMAX  -  maximum  velocity  occurring  at  an  outflow  boundary 

UGLIM  -  velocity  returned  by  subroutine  GLIMM 

UL  -  left  side  velocity  for  RP 

UR  -  right  side  velocity  for  RP 

USTAR  -  velocity  in  'starred'  state  of  RP  solution  (see  Fig.  3) 

WDP  -  value  returned  by  random  number  generator  subprogram 

WL  -  left  shock  wave  velocity 

WR  -  right  shock  wave  velocity 

WNORM  -  variable  used  in  random  number  generator  subprogram 

X  -  space  dimension 

XCNTCT  -  location  of  contact  surface 

XI, XII  -  random  numbers  scaled  to  grid  cell 
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XREF     -  reference  length 

Y        -  argument  used  in  function  subprogram  PHI  equal  to  PSTAR 

Z        -  argument  used  in  function  subprogram  PHI  equal  to  either 

PL  or  PR 
ZETA     -  dummy  variable  (for  initialization  purposes  in  random 

number  generator) 


B. 5.   List  of  Subroutines,  Function  Subprograms 


B.5.1.   Subroutines 


INIT1    -  prescribes  initial  data  corresponding  to  SWL=1,  SWR= 1 ; 
e.g.,  shock-tube  problem 

INIT2L  -  prescribes  initial  data  corresponding  to  SWL=2,  SWR=1 

INIT2R  -  prescribes  initial  data  corresponding  to  SWL=1,  SWR=2 

INIT3L  -  prescribes  initial  data  corresponding  to  SWL=3,  SWR=1 

INIT3R  -  prescribes  initial  data  corresponding  to  SWL=1,  SWR=3 

PL0Ti,2  -  graphics  subroutines 

GLIMM   -  solves  the  Riemann  problem,  samples  the  solution  and 
returns  values  for  flow  parameters 


GE 


modular  user  supplied  subroutine  to  simulate  wave  diagram 
of  General  Electric  Wave  Engine 


DETON   -  modular  user  supplied  subroutine  to  simulate  evacuation  of 
detonation  chamber 

SPCTRA  -  modular  user  supplied  subroutine  to  simulate  wave  diagram 
of  Spectra  Technology's  Pressure  Exchanger 

BCL1     -  prescribes  boundary  conditions  (BC's)  corresponding  to 
SWL=1,  i.e.,  solid  wall  on  left  side 

BCL2    -  prescribes  BC's  corresponding  to  SWL=2,  i.e.,  outflow  at 
constant  static  pressure  on  left  side 

BCL3    -  prescribes  BC's  corresponding  to  SWL=3,  i.e.,  'piston' 
inflow  on  left  side 

BCL4    -  prescribes  BC's  corresponding  to  SWL=4,  i.e.,  isentropic 
inflow  from  reservoir  on  left  side 

BCL5    -  prescribes  BC's  corresponding  to  SWL=5,  i.e.,  wave  'tuning' 
on  left  side 
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8CR1,  BCR2,  BCR3,  BCR4 ,  BCR5  -  prescribe  BC's  corresponding  to 
SWR=1 ,2,3,4,5  respectively  on  right  side 


B.5.2.   Function  Subprograms 

PHI(y,z)   -  required  in  iteration  procedure  for  solution  of  RP 

PHIl(PB)   -  describes  shock  transition  function,   ^a(Pb)  »  f°r 
two  states  a  and  b  connected  by  a  shock  wave  (see 
Ref.  6,  Ch.  Ill) 

PSI(PB)    -  describes  rarefaction  transition  function,   ^a(Pb)   » 
for  two  states  a  and  b  connected  by  a  rarefaction  wave 
(see  Ref.  6,  Ch.  Ill) 

WDP(II)   -  generates  a  random  number  in  a  van  der  Corput  sequence 

each  time  it  is  invoked.   Note  that  it  needs  to  be  called 
once  from  outside  the  main  loop  by  specifying  an  argument 
11=1  to  initialize  IDIGT  and  WNORM,  returning  a  value 
of  0  for  the  dummy  variable  ZETA,  and  then  a  second  time 
from  within  the  main  loop  with  an  argument  11=0  to  return 
a  value  which  is  the  random  number. 
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